{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3561e8aa",
   "metadata": {},
   "source": [
    "# Gauss curvature approximation on embedded surfaces\n",
    "\n",
    "In the [previous notebook](distributional_gauss_curvature.ipynb), we presented and discussed the distributional Gauss curvature and its lifting on Riemannian manifolds. However, we can also apply the method of high-order Gauss curvature on a surface $\\mathcal{S}$ embedded in $\\mathbb{R}^3$ with only a few adaptions of the method. Note that the theory can partly be extended from Riemannian manifolds to surfaces. \n",
    "\n",
    "## Setting\n",
    "The surface can be expressed by its embedding $\\Phi$, which is a polynomial. If piecewise linear, we obtain a linear/affine approximation of the exact surface. If it is of higher order, the discrete geometry gets curved for better approximation. From $\\Phi$ we can compute its gradient $F=\\nabla \\Phi\\in \\mathbb{R}^{3\\times 2}$ and the metric tensor induced by the embedding via $g= F^TF\\in \\mathbb{R}^{2\\times 2}$. Note that the polynomial degree of $g$ is $2(k-1)$ if $k$ is the order of the embedding $\\Phi$. Although the theory is available for such a $g$, we won't use a discrete lifted Gauss curvature of order $2(k-1)+1$ to fit the theory, but the more natural choice of order $k$. Furthermore, to make use of the theory that for a discrete metric stemming from the canonical Regge interpolant, we would need a commuting property of the form $\\mathcal{I}^{\\mathrm{Reg}^{2(k-1)}}(g)= (\\nabla \\mathcal{I}^{\\mathrm{Lag}^k}\\Phi)^T(\\nabla \\mathcal{I}^{\\mathrm{Lag}^k}\\Phi)$, which in general does not hold.\n",
    "\n",
    "\n",
    "Let $\\mathcal{S}$ be a surface with triangulation $\\mathcal{T}_h$.  The finite element method to approximate the Gauss curvature reads: Find $K_h \\in V_h^k$ (the Lagrange finite element space of degree $k$) such that for all $\\varphi \\in V_h^k$, \n",
    "$$\n",
    "\\int_{\\mathcal{T}_h} K_h \\varphi \\,da = \\sum_{T \\in \\mathcal{T}_h} \\int_T K|_T \\varphi\\, da + \\sum_{E \\in \\mathcal{E}_h} \\int_E [\\![ \\kappa_g ]\\!] \\varphi\\, d\\ell+ \\sum_{V \\in \\mathcal{V}_h} \\Theta_V \\varphi(V).\n",
    "$$\n",
    "Here $K|_T$ is the Gauss curvature in the interior of a triangle $T$, $[\\![\\kappa_g]\\!]$ is the jump of the geodesic curvature $\\kappa_g = \\mu \\cdot \\nabla_t t$ across the edge $E$, and \n",
    "$$\n",
    "\\Theta_V = 2\\pi-\\sum_{T:V\\subset T}\\sphericalangle_V^T\n",
    "$$\n",
    "denotes the angle defect at the vertex $V$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7a1fc3c",
   "metadata": {},
   "source": [
    "## Numerical example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c455fe7",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from netgen.occ import *\n",
    "from ngsolve.webgui import Draw\n",
    "import random as random\n",
    "from ngsolve.krylovspace import CGSolver"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4ac5eb4",
   "metadata": {},
   "source": [
    "We define a function that computes the $H^{-1}$-norm of a functional $f$, using the fact that $||f||_{H^{-1}}$ is equivalent to $\\|u\\|_{H^1}$ if $-\\Delta u = f$.  We will use this function later to compute the error in the discrete Gauss curvature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19e2ee13",
   "metadata": {},
   "outputs": [],
   "source": [
    "# H^-1 norm\n",
    "def ComputeHm1Norm(rhs, order):\n",
    "    fesH = H1(mesh, order=order)\n",
    "    u, v = fesH.TnT()\n",
    "\n",
    "    a = BilinearForm(\n",
    "        (Grad(u).Trace() * Grad(v).Trace() + u * v) * ds,\n",
    "        symmetric=True,\n",
    "        symmetric_storage=True,\n",
    "        condense=True,\n",
    "    )\n",
    "    f = LinearForm(rhs * v * ds).Assemble()\n",
    "\n",
    "    apre = Preconditioner(a, \"bddc\")\n",
    "    a.Assemble()\n",
    "    invS = CGSolver(a.mat, apre.mat, printrates=\"\\r\", maxiter=400)\n",
    "    ext = IdentityMatrix() + a.harmonic_extension\n",
    "    inv = a.inner_solve + ext @ invS @ ext.T\n",
    "\n",
    "    gfu = GridFunction(fesH)\n",
    "    gfu.vec.data = inv * f.vec\n",
    "\n",
    "    err = sqrt(Integrate(gfu**2 + Grad(gfu) ** 2, mesh, BND))\n",
    "    return err"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9927b3da",
   "metadata": {},
   "source": [
    "Next we define a function that computes $K_h = \\sum_i u_i \\varphi_i$ by solving the linear system $Mu=f$, where $M_{ij} = \\int_{\\mathcal{T}_h} \\varphi_i \\varphi_j da$ and \n",
    "$$\n",
    "f_i = \\sum_{T \\in \\mathcal{T}_h} \\int_T K|_T \\varphi_i\\, da+ \\sum_{E \\in \\mathcal{E}_h} \\int_E [\\![ \\kappa_g ]\\!] \\varphi_i \\,d\\ell  +\\sum_{V \\in \\mathcal{V}_h} \\Theta_V \\varphi_i(V).\n",
    "$$\n",
    "After computing $K_h$, the function outputs the errors $||K_h-K||_{H^{-1}}$ and $||K_h-K||_{L^2}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5d02c26",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ComputeGaussCurvature(mesh, order, Kex):\n",
    "    fes = H1(mesh, order=order)\n",
    "    u, v = fes.TnT()\n",
    "\n",
    "    # for angle deficit\n",
    "    bbnd_tang = specialcf.VertexTangentialVectors(3)\n",
    "    bbnd_tang1 = bbnd_tang[:, 0]\n",
    "    bbnd_tang2 = bbnd_tang[:, 1]\n",
    "    # for geodesic curvature\n",
    "    mu = Cross(specialcf.normal(3), specialcf.tangential(3))\n",
    "    edge_curve = specialcf.EdgeCurvature(3)  # nabla_t t\n",
    "\n",
    "    # for elementwise Gauss curvature\n",
    "    def GaussCurvature():\n",
    "        nsurf = specialcf.normal(3)\n",
    "        return Cof(Grad(nsurf)) * nsurf * nsurf\n",
    "\n",
    "    # distributional Gauss curvature\n",
    "    f = LinearForm(fes)\n",
    "    # elementwise Gauss curvature\n",
    "    f += GaussCurvature() * v * ds\n",
    "    # jump of geodesic curvature\n",
    "    f += -edge_curve * mu * v * ds(element_boundary=True)\n",
    "    # one part of angle defect\n",
    "    f += -v * acos(bbnd_tang1 * bbnd_tang2) * ds(element_vb=BBND)\n",
    "\n",
    "    # mass matrix to compute discrete L2 Riesz representative\n",
    "    M = BilinearForm(u * v * ds, symmetric=True, symmetric_storage=True, condense=True)\n",
    "\n",
    "    gf_K = GridFunction(fes)\n",
    "\n",
    "    f.Assemble()\n",
    "    # second part of angle deficit (closed surface, no boundary)\n",
    "    for i in range(mesh.nv):\n",
    "        f.vec[i] += 2 * pi\n",
    "\n",
    "    Mpre = Preconditioner(M, \"bddc\")\n",
    "    M.Assemble()\n",
    "    invS = CGSolver(M.mat, Mpre.mat, printrates=\"\\r\", maxiter=400)\n",
    "    ext = IdentityMatrix() + M.harmonic_extension\n",
    "    inv = M.inner_solve + ext @ invS @ ext.T\n",
    "\n",
    "    gf_K.vec.data = inv * f.vec\n",
    "\n",
    "    l2_err = sqrt(Integrate((gf_K - Kex) ** 2, mesh, BND))\n",
    "    hm1_err = ComputeHm1Norm(gf_K - Kex, order=order + 2)\n",
    "    print(\n",
    "        f\" Check Gauss-Bonnet theorem: int gf_K = {Integrate(gf_K * ds(bonus_intorder=5), mesh)}\"\n",
    "        f\" = 4*pi = {4 * pi}\"\n",
    "    )\n",
    "\n",
    "    # uncomment to draw Gauss curvature\n",
    "    Draw(gf_K, mesh, \"K\")\n",
    "\n",
    "    return l2_err, hm1_err, fes.ndof"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9cc2fdb",
   "metadata": {},
   "source": [
    "Now, we are ready to test the method. We triangulate an ellipsoid (using curved triangles of a given order), compute the Gaussian curvature, measure the error, and repeat several refinements. We also check how well the lifted Gauss curvature fulfils the Gauss-Bonnet theorem $\\int_{\\mathcal{T}_h}K_h\\,da = 2\\pi\\chi(\\mathcal{S}) = 4\\pi$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0e1f796",
   "metadata": {},
   "outputs": [],
   "source": [
    "order = 1\n",
    "# radius of sphere/ellipsoid\n",
    "R = 3\n",
    "\n",
    "# use sphere or ellipsoid surface\n",
    "if False:  # sphere\n",
    "    geo = Sphere((0, 0, 0), R).faces[0]\n",
    "    # exact Gauss curvature of sphere\n",
    "    Kex = 1 / R**2\n",
    "else:  # ellipsoid\n",
    "    a = R\n",
    "    b = R\n",
    "    c = 3 / 4 * R\n",
    "    geo = Ellipsoid(Axes((0, 0, 0), X, Y), a, b, c).faces[0]\n",
    "    # exact Gauss curvature of ellipsoid\n",
    "    Kex = 1 / (a**2 * b**2 * c**2 * (x**2 / a**4 + y**2 / b**4 + z**2 / c**4) ** 2)\n",
    "\n",
    "err_l2 = []\n",
    "err_hm1 = []\n",
    "ndof = []\n",
    "\n",
    "with TaskManager():\n",
    "    for i in range(4 + (order == 1)):\n",
    "        mesh = Mesh(OCCGeometry(geo).GenerateMesh(maxh=0.5**i)).Curve(order)\n",
    "\n",
    "        errl, errm1, dof = ComputeGaussCurvature(mesh, order, Kex=Kex)\n",
    "        err_l2.append(errl)\n",
    "        err_hm1.append(errm1)\n",
    "        ndof.append(dof)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f7fc6ec",
   "metadata": {},
   "source": [
    "Finally, we plot the errors against the number of degrees of freedom."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51a457f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.plot(ndof, err_hm1, \"-o\", label=r\"$\\|K_h-K_{\\mathrm{ex}}\\|_{H^{-1}}$\")\n",
    "plt.plot(ndof, err_l2, \"-x\", label=r\"$\\|K_h-K_{\\mathrm{ex}}\\|_{L^2}$\")\n",
    "plt.plot(\n",
    "    ndof,\n",
    "    [order**order * dof ** (-(order - 1) / 2) for dof in ndof],\n",
    "    \"-.\",\n",
    "    color=\"k\",\n",
    "    label=f\"$O(h^{order-1})$\",\n",
    ")\n",
    "plt.plot(\n",
    "    ndof,\n",
    "    [order**order * dof ** (-(order) / 2) for dof in ndof],\n",
    "    \"-\",\n",
    "    color=\"k\",\n",
    "    label=f\"$O(h^{order})$\",\n",
    ")\n",
    "plt.plot(\n",
    "    ndof,\n",
    "    [order**order * dof ** (-(order + 1) / 2) for dof in ndof],\n",
    "    \"--\",\n",
    "    color=\"k\",\n",
    "    label=f\"$O(h^{order+1})$\",\n",
    ")\n",
    "plt.yscale(\"log\")\n",
    "plt.xscale(\"log\")\n",
    "plt.legend()\n",
    "plt.title(f\"order = {order}\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0678ecce",
   "metadata": {},
   "source": [
    "Observations: Except for $k=2$, we observe the convergence rates\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\|K_h-K\\|_{L^2}\\le C h^{k-1},\\qquad \\|K_h-K\\|_{H^{-1}}\\le C h^{k}.\n",
    "\\end{align*}\n",
    "$$\n",
    "One additional order of convergence is observed for $k=2$ of curving the surface quadratically. This could be a lucky case where the theory of improved convergence rates applies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da0c7101",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}