{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3561e8aa",
   "metadata": {},
   "source": [
    "# Distributional scalar curvature\n",
    "After the [distributional Gauss curvature](distributional_gauss_curvature.ipynb), we consider the distributional scalar curvature for arbitrarily dimensional Riemannian manifolds in this notebook. In two dimensions, the scalar curvature and the Gauss curvature coincide up to a factor $2$. We will see that the angle defect directly extends to a higher dimension and that the geodesic curvature becomes the mean curvature. Further, the integral representation of the distributional scalar curvature splits into two parts, where the second one is zero in the two-dimensional case and can be identified with the distributional Einstein tensor discussed in the [next notebook](distributional_einstein_tensor.ipynb). For a detailed derivation and analysis, we refer to [Gawlik, Neunteufel. *Finite element approximation of scalar curvature in arbitrary dimension*, *Math. Comp.* (2024).](https://doi.org/10.1090/mcom/4038).\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "## Scalar curvature as a nonlinear distribution\n",
    "The scalar curvature $S$ is the double contraction of the Riemann curvature tensor $\\mathfrak{R}$ or equivalently the contraction of the Ricci curvature tensor $\\mathrm{Ric}$\n",
    "\\begin{align*}\n",
    "S = g_{ij}g_{kl}\\mathfrak{R}_{iklj}=g_{ij}\\mathrm{Ric}_{ij}.\n",
    "\\end{align*}\n",
    "With the same arguments as for the Gauss curvature, the scalar curvature is a nonlinear (quasi-linear) second-order differential operator of the metric. \n",
    "\n",
    "The distributional (densitized) scalar curvature acting on a sufficiently smooth scalar-valued test function $\\varphi$ is defined as\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\widetilde{S\\omega}(\\varphi) :=\\sum_{T\\in\\mathcal{T}}\\int_TS|_T\\,\\varphi\\,\\omega_T+2\\sum_{F\\in\\mathcal{F}}\\int_F[\\![H]\\!]\\,\\varphi\\,\\omega_F+2\\sum_{E\\in\\mathcal{E}}\\int_E\\Theta_E\\,\\varphi\\,\\omega_E,\n",
    "\\end{align*}\n",
    "$$\n",
    "where $\\omega_T$, $\\omega_F$, and $\\omega_E$ are the respective volume forms on elements $T$, co-dimension 1 facets $F$ and co-dimension two boundaries $E$ (vertices in 2D, edges in 3D, etc.). $H$ denotes the mean curvature defined as the trace of the second fundamental form $I\\!I$\n",
    "$$\n",
    "\\begin{align*}\n",
    "H:= \\mathrm{tr}(I\\!I)= g_{ij}I\\!I^{ij},\\qquad I\\!I(X,Y) = g(\\nabla_Xn,Y)=-g(\\nabla_XY,n),\\quad X,Y\\in \\mathfrak{X}(F),\n",
    "\\end{align*}\n",
    "$$\n",
    "where $n$ is the normal vector at the facet (caution: there exist again different sign conventions), and $\\Theta_E$ is the angle defect. Note that the concept of the angle defect directly extends from 2D to any dimension, as we always measure the angle in the two-dimensional plane perpendicular to the co-dimension 2 boundary. Thus, exactly the same formula applies. \n",
    "\n",
    "Although there is no Gauss-Bonnet theorem in arbitrary dimension, we can heuristically motivate the jump of the mean curvature as the appropriate facet term. At a point of a facet $p\\in F$, the scalar curvature can be expressed as the sum of sectional curvatures of $N(N-1)/2$ tangent planes, $N$ denoting the dimension, that are pairwise $g$-perpendicular at $p$. $(N-1)(N-2)/2$ of those planes are tangent to $F$ and $N-1$ ones are $g$-perpendicular to $F$ at $p$. Because of the tangential-tangential continuity of the Regge metric, the sectional curvature of the planes tangent to $F$ are non-singular and do not yield a contribution to the facet term. The other $N-1$ planes, however, have jumping curvature. On these two-dimensional planes, we can apply the Gauss-Bonnet theorem. Summing up those ``geodesic curvature jumps'' yields the jump of the mean curvature of $F$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91fc6eab",
   "metadata": {},
   "source": [
    "## Integral representation and convergence\n",
    "\n",
    "Taking the variation of the densitized distributional scalar curvature in the metric in direction $\\sigma$ yields\n",
    "$$\n",
    "\\begin{align*}\n",
    "D_g\\widetilde{S\\omega}(v)[\\sigma] = b_h(g;\\sigma,v)-a_h(g;\\sigma,v)\n",
    "\\end{align*}\n",
    "$$\n",
    "with\n",
    "$$\n",
    "\\begin{align*}\n",
    "b_h(g;\\sigma,v)&= \\sum_T\\int_T\\langle\\mathbb{S}\\sigma,\\nabla^2v\\rangle\\,\\omega_T-\\sum_{F}\\int_F\\mathbb{S}\\sigma(n,n)[\\![\\nabla_n v]\\!]_F\\,\\omega_F,\\\\\n",
    "a_h(g;\\sigma,v)&= \\sum_T\\int_T\\langle G,\\sigma\\rangle\\,v\\,\\omega_T+\\sum_F\\int_F\\langle [\\![\\bar{I\\!I}]\\!]_F,\\sigma|_F\\rangle\\,v\\omega_F-\\sum_{E}\\int_E\\langle\\Theta_Eg|_E,\\sigma|_E\\rangle\\,v\\,\\omega_E.\n",
    "\\end{align*}\n",
    "$$\n",
    "Here, all differential operators are covariant with respect to the metric $g$. $\\bar{I\\!I}=I\\!I-H\\,g|_F$ denotes the trace-reversed second fundamental form and $G=\\mathrm{Ric}-1/2S\\,g$ is the Einstein tensor. Note that in two dimensions, $a(g;\\sigma,v)$ becomes identically zero, and $b_h(g;\\sigma,v)$ coincides with the expression of the integral representation of the distributional Gauss curvature. As we will discuss in the [following notebook](distributional_einstein_tensor.ipynb), $a_h(g;\\sigma,v)$ can be interpreted as the distributional Einstein tensor when we set $v=1$. $b_h(g;\\sigma,v)$ can also be interpreted as extending the (covariant) Hellan-Herrmann-Johnson method to arbitrary dimensions. Its numerical analysis follows the same lines as for the Gauss curvature. In three and higher dimensions, however, the bilinear form $a_h(g;\\sigma,v)$ needs to be estimated too. \n",
    "\n",
    "We obtain the following convergence rates in the $H^{-2}$ norm for $k\\ge 0$ and $g_h\\in \\mathrm{Reg}^k$\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\|\\widetilde{S\\omega}(g_h)-(S\\omega)(g)\\|_{H^{-2}}\\le C\\begin{cases}\n",
    "h^{k+1}, & N=2\\\\\n",
    "h^{\\min\\{k+1,2k\\}}, &N \\ge 3\n",
    "\\end{cases}.\n",
    "\\end{align*}\n",
    "$$\n",
    "We proved that in 2D, the angle defect for approximating the scalar (and therefore Gauss) curvature always converges in the $H^{-2}$ norm. Before, only convergence in the sense of measures was known. We cannot expect convergence in the lowest-order case, $k=0$, of a piecewise constant Regge metric for three or higher dimensions. When going to linear order, the convergence rate jumps to be quadratic. Then, the dependency on the polynomial approximation order relates to the convergence rate as expected. The co-dimension 2 term of $a_h(g;\\sigma,v)$ turns out to be the issue for not obtaining convergence in the lowest-order case. Unfortunately, numerical investigations show that this estimate is sharp. Nevertheless, in several examples, we still observed linear convergence as the absolute error of this term seems to be small, and the error of the other terms dominates such that the stagnation of error is hard to spot. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43e1528e",
   "metadata": {},
   "source": [
    "## Numerical example\n",
    "In the following example, we compute the distributional scalar curvature $\\widetilde{S\\omega}$ on a sequence of meshes for a three-dimensional Riemannian manifold and compute the $H^{-2}$ error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "daf1e654",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve.meshes import MakeStructured3DMesh\n",
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "import random as random\n",
    "from ngsolve.krylovspace import CGSolver\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "def mapping(x, y, z):\n",
    "    return (2 * x - 1, 2 * y - 1, 2 * z - 1)\n",
    "\n",
    "\n",
    "mesh = MakeStructured3DMesh(False, nx=4, ny=4, nz=4, mapping=mapping)\n",
    "\n",
    "# Define an embedded \"3D-surface\" in 4D via a graph function\n",
    "peak = (\n",
    "    0.5 * x**2 - 1 / 12 * x**4 + 0.5 * y**2 - 1 / 12 * y**4 + 0.5 * z**2 - 1 / 12 * z**4\n",
    ")\n",
    "\n",
    "# compute the exact induced metric tensor G\n",
    "F = CF(\n",
    "    (1, 0, 0, 0, 1, 0, 0, 0, 1, peak.Diff(x), peak.Diff(y), peak.Diff(z)), dims=(4, 3)\n",
    ")\n",
    "G_ex = F.trans * F\n",
    "\n",
    "\n",
    "# exact scalar curvature\n",
    "def q(x):\n",
    "    return x**2 * (x**2 - 3) ** 2\n",
    "\n",
    "\n",
    "# exact scalar curvature\n",
    "S_ex = (\n",
    "    18\n",
    "    * (\n",
    "        (1 - x**2) * (1 - y**2) * (9 + q(z))\n",
    "        + (1 - y**2) * (1 - z**2) * (9 + q(x))\n",
    "        + (1 - z**2) * (1 - x**2) * (9 + q(y))\n",
    "    )\n",
    ") / ((9 + q(x) + q(y) + q(z)) ** 2)\n",
    "\n",
    "Draw(S_ex, mesh);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d35a8ac3",
   "metadata": {},
   "source": [
    "We compute the $H^{-2}$-norm of the error $\\|\\widetilde{S\\omega}-S\\omega\\|_{H^{-2}}$ by solving the following fourth-order problem $\\mathrm{div}(\\mathrm{div}(\\nabla^2 u))= \\widetilde{S\\omega}-S\\omega$. Then $\\|u\\|_{H^2}\\approx \\|\\widetilde{S\\omega}-S\\omega\\|_{H^{-2}}$. We use the Euclidean Hellan-Herrmann-Johnson (HHJ) method to solve the fourth-order problem as a mixed saddle point problem. Therein, $\\sigma=\\nabla^2u$ is introduced as an additional unknown. We use hybridization techniques to regain a minimization problem. We use a BDDC preconditioner in combination with static condensation and a CGSolver to solve the problem iteratively. To mitigate approximation errors that spoil the convergence, we use two polynomial degrees more than for $g_h$ to compute $u_h$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1e60dfd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute H^-2 error via HHJ method\n",
    "def ComputeHm2Norm(func_vol, func_bnd, func_bbnd, mesh, order):\n",
    "    fesH = H1(mesh, order=order, dirichlet=\".*\")\n",
    "    fesS = Discontinuous(HDivDiv(mesh, order=order - 1))\n",
    "    fesHB = NormalFacetFESpace(mesh, order=order - 1, dirichlet=\".*\")\n",
    "    fes = fesH * fesS * fesHB\n",
    "    (u, sigma, uh), (v, tau, vh) = fes.TnT()\n",
    "\n",
    "    n = specialcf.normal(mesh.dim)\n",
    "\n",
    "    a = BilinearForm(fes, symmetric=True, symmetric_storage=True, condense=True)\n",
    "    a += (\n",
    "        InnerProduct(sigma, tau)\n",
    "        - InnerProduct(tau, u.Operator(\"hesse\"))\n",
    "        - InnerProduct(sigma, v.Operator(\"hesse\"))\n",
    "    ) * dx\n",
    "    a += ((Grad(u) - uh) * n * tau[n, n] + (Grad(v) - vh) * n * sigma[n, n]) * dx(\n",
    "        element_boundary=True\n",
    "    )\n",
    "\n",
    "    f = LinearForm(fes)\n",
    "    f += (\n",
    "        func_vol * v * dx\n",
    "        + func_bnd * v * dx(element_boundary=True)\n",
    "        + func_bbnd * v * dx(element_vb=BBND)\n",
    "    )\n",
    "    f.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",
    "    gf_error_lift = GridFunction(fes)\n",
    "    gf_error_lift.vec.data = inv * f.vec\n",
    "\n",
    "    gf_u, gf_sigma, _ = gf_error_lift.components\n",
    "    # \"sigma = hesse u\"\n",
    "    err = sqrt(\n",
    "        Integrate(gf_u**2 + Grad(gf_u) ** 2 + InnerProduct(gf_sigma, gf_sigma), mesh)\n",
    "    )\n",
    "    return err"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74acccd7",
   "metadata": {},
   "source": [
    "The following function computes the $H^2$-error of the distributional scalar curvature for a given metric and polynomial order. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5f0e469",
   "metadata": {},
   "outputs": [],
   "source": [
    "# main function to compute the scalar curvature\n",
    "def ComputeScalarCurvature(mesh, G, order):\n",
    "    # utility vectors for the codimension 2 boundary\n",
    "    bbnd_tang = specialcf.EdgeFaceTangentialVectors(3)\n",
    "    tE = specialcf.tangential(mesh.dim, True)\n",
    "    tef1 = bbnd_tang[:, 0]\n",
    "    tef2 = bbnd_tang[:, 1]\n",
    "    n1 = Cross(tef1, tE)\n",
    "    n2 = Cross(tef2, tE)\n",
    "    # facet normal vector\n",
    "    n = specialcf.normal(mesh.dim)\n",
    "\n",
    "    # Regge finite element space\n",
    "    fesCC = HCurlCurl(mesh, order=order)\n",
    "    gfG = GridFunction(fesCC)\n",
    "    gfG.Set(G)\n",
    "\n",
    "    # volume forms\n",
    "    vol_bbnd = sqrt(gfG * tE * tE)\n",
    "    vol_bnd = sqrt(Cof(gfG) * n * n)\n",
    "    vol = sqrt(Det(gfG))\n",
    "\n",
    "    # angle deficit\n",
    "    n1g = 1 / sqrt(Inv(gfG) * n1 * n1) * Inv(gfG) * n1\n",
    "    n2g = 1 / sqrt(Inv(gfG) * n2 * n2) * Inv(gfG) * n2\n",
    "    func_bbnd = 2 * (acos(n1 * n2) - acos(gfG * n1g * n2g)) * vol_bbnd\n",
    "\n",
    "    # jump of mean curvature\n",
    "    ng = 1 / sqrt(Inv(gfG) * n * n) * Inv(gfG) * n\n",
    "    SFF = -gfG.Operator(\"christoffel\") * ng\n",
    "    P = Inv(gfG) - OuterProduct(ng, ng)\n",
    "    H = InnerProduct(gfG, P * SFF * P)\n",
    "    func_bnd = 2 * H * vol_bnd\n",
    "    # scalar curvature minus exact one\n",
    "    func_vol = gfG.Operator(\"scalar\") * vol - S_ex * sqrt(Det(G))\n",
    "\n",
    "    err_hm2 = ComputeHm2Norm(func_vol, func_bnd, func_bbnd, mesh, order=order + 2)\n",
    "\n",
    "    return (err_hm2, fesCC.ndof)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb991c32",
   "metadata": {},
   "source": [
    "To avoid possible super-convergence due to symmetries in the mesh, we perturb the internal vertices by a random noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9643f6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute the diameter of the mesh as maximal mesh size h\n",
    "def CompDiam(mesh):\n",
    "    lengths = np.zeros(mesh.nedge)\n",
    "    for edge in mesh.edges:\n",
    "        pnts = []\n",
    "        for vert in edge.vertices:\n",
    "            pnt = mesh[vert].point\n",
    "            pnts.append((pnt[0], pnt[1], pnt[2]))\n",
    "        lengths[edge.nr] = sqrt(\n",
    "            (pnts[0][0] - pnts[1][0]) ** 2\n",
    "            + (pnts[0][1] - pnts[1][1]) ** 2\n",
    "            + (pnts[0][2] - pnts[1][2]) ** 2\n",
    "        )\n",
    "\n",
    "    print(f\"max = {max(lengths)}, min = {min(lengths)}, ne = {mesh.ne}\")\n",
    "    return max(lengths)\n",
    "\n",
    "\n",
    "# generate and randomly perturb it, check if the resulting mesh is valid\n",
    "def GenerateRandomMesh(k, par):\n",
    "    for i in range(20):\n",
    "        mesh = MakeStructured3DMesh(False, nx=2**k, ny=2**k, nz=2**k, mapping=mapping)\n",
    "        maxh = CompDiam(mesh)\n",
    "        for pnts in mesh.ngmesh.Points():\n",
    "            px, py, pz = pnts[0], pnts[1], pnts[2]\n",
    "            if (\n",
    "                abs(px + 1) > 1e-8\n",
    "                and abs(px - 1) > 1e-8\n",
    "                and abs(py + 1) > 1e-8\n",
    "                and abs(py - 1) > 1e-8\n",
    "                and abs(pz + 1) > 1e-8\n",
    "                and abs(pz - 1) > 1e-8\n",
    "            ):\n",
    "                px += random.uniform(-maxh / 2**par / sqrt(2), maxh / 2**par / sqrt(2))\n",
    "                py += random.uniform(-maxh / 2**par / sqrt(2), maxh / 2**par / sqrt(2))\n",
    "                pz += random.uniform(-maxh / 2**par / sqrt(2), maxh / 2**par / sqrt(2))\n",
    "            pnts[0] = px\n",
    "            pnts[1] = py\n",
    "            pnts[2] = pz\n",
    "        F = specialcf.JacobianMatrix(3)\n",
    "        # check if the mesh is valid\n",
    "        if Integrate(IfPos(Det(F), 0, -1), mesh) >= 0:\n",
    "            break\n",
    "    else:\n",
    "        raise Exception(\"Could not generate a valid mesh\")\n",
    "\n",
    "    return mesh"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df00a443",
   "metadata": {},
   "source": [
    "Compute the scalar curvature approximation on a sequence of meshes. You can play around with the polynomial order. The $H^{-2}$ error and the number of degrees of freedom are stored."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0080abbf",
   "metadata": {},
   "outputs": [],
   "source": [
    "err_hm2 = []\n",
    "ndof = []\n",
    "\n",
    "# order of metric approximation >= 0\n",
    "order = 1\n",
    "\n",
    "with TaskManager():\n",
    "    for i in range(0, 4):\n",
    "        mesh = GenerateRandomMesh(k=i, par=3)\n",
    "        errm2, dof = ComputeScalarCurvature(mesh, G=G_ex, order=order)\n",
    "        err_hm2.append(errm2)\n",
    "        ndof.append(dof)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "242ee313",
   "metadata": {},
   "source": [
    "Plot the results together with reference lines for linear, quadratic, and cubic convergence order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0282a527",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.plot(\n",
    "    ndof,\n",
    "    err_hm2,\n",
    "    \"-x\",\n",
    "    label=r\"$\\|\\widetilde{S\\omega}-S_{\\mathrm{ex}}\\omega_{\\mathrm{ex}}\\|_{H^{-2}}$\",\n",
    ")\n",
    "plt.plot(ndof, [2 * dof ** (-1 / 3) for dof in ndof], \"-\", color=\"k\", label=\"$O(h)$\")\n",
    "plt.plot(ndof, [5 * dof ** (-2 / 3) for dof in ndof], \"--\", color=\"k\", label=\"$O(h^2)$\")\n",
    "plt.plot(\n",
    "    ndof, [50 * dof ** (-3 / 3) for dof in ndof], \"-.\", color=\"k\", label=\"$O(h^3)$\"\n",
    ")\n",
    "\n",
    "\n",
    "plt.yscale(\"log\")\n",
    "plt.xscale(\"log\")\n",
    "plt.legend()\n",
    "plt.title(f\"order = {order}\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea93ff47",
   "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
}