{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3561e8aa",
   "metadata": {},
   "source": [
    "# Distributional Gauss curvature (analysis)\n",
    "After we defined the distributional Gauss curvature and its lifting and observed convergence rates numerically in the [previous notebook](distributional_gauss_curvature.ipynb), we focus in this notebook on the numerical analysis and presenting ideas on how to prove the observed convergence rates.\n",
    "\n",
    "The notebook is based on the following works: [Gawlik. *High-Order Approximation of Gaussian Curvature with Regge Finite Elements*, *SIAM Journal on Numerical Analysis* (2020).](https://doi.org/10.1137/19M1255549), [Berchenko-Kogan, Gawlik. *Finite Element Approximation of the Levi-Civita Connection and Its Curvature in Two Dimensions*, *Found Comput Math* (2022)](https://doi.org/10.1007/s10208-022-09597-1), [Gopalakrishnan, Neunteufel, Schöberl, Wardetzky. *Analysis of curvature approximations via covariant curl and incompatibility for Regge metrics*, *The SMAI Journal of computational mathematics* (2023).](https://doi.org/10.5802/smai-jcm.98), and [Gopalakrishnan, Neunteufel, Schöberl, Wardetzky. *On the improved convergence of lifted distributional Gauss curvature from Regge elements*, *Results in Applied Mathematics* (2024).](https://doi.org/10.1016/j.rinam.2024.100511).\n",
    "\n",
    "The distributional Gauss curvature\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\widetilde{K\\,\\omega}(\\varphi) = \\sum_{T\\in\\mathcal{T}}\\int_TK|_T\\,\\varphi\\,\\omega|_T+\\sum_{E\\in\\mathcal{E}}\\int_E[\\![\\kappa_g]\\!]\\,\\varphi\\,\\omega_E+\\sum_{V\\in\\mathcal{V}}\\sphericalangle_V\\,\\varphi(V)\n",
    "\\end{align*}\n",
    "$$\n",
    "is highly nonlinear in the metric, making numerical analysis unfeasible in its current form. Fortunately, the distributional Gauss curvature entails an integral representation generalizing the Hellan-Herrmann-Johnson method from the Euclidean to the covariant setting."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca0929bc",
   "metadata": {},
   "source": [
    "## Integral representation of distributional Gauss curvature\n",
    "\n",
    "Christiansen showed in [Christiansen. *On the linearization of Regge calculus*. *Numer. Math.* (2011).](https://doi.org/10.1007/s00211-011-0394-z) that the linearization of the Einstein-Hilbert integral around the Euclidean metric is the incompatibility operator in 3D. Gawlik investigated in [Gawlik. *High-Order Approximation of Gaussian Curvature with Regge Finite Elements*, *SIAM Journal on Numerical Analysis* (2020).](https://doi.org/10.1137/19M1255549) the variation of the densitized Gauss curvature around a given general metric. There holds for a smooth metric the following Gateaux derivative in a direction $\\sigma\\in C^{\\infty}(\\Omega,\\mathbb{R}^{2\\times 2}_{\\mathrm{sym}})$\n",
    "\\begin{align*}\n",
    "D_g(K(g)\\omega(g))[\\sigma]:= \\lim\\limits_{t\\to 0}\\frac{K(g+t\\sigma)\\omega(g+t\\sigma)-K(g)\\omega(g)}{t} = \\frac{1}{2}\\mathrm{div}_g\\mathrm{div}_g(\\mathbb{S}_{g}\\sigma)\\omega(g),\n",
    "\\end{align*}\n",
    "where $\\mathbb{S}_g\\sigma = \\sigma -\\mathrm{tr}_g(\\sigma)g$ and subscript $g$ indicates that [covariant differential operators](covariant_derivatives.ipynb) are considered.\n",
    "There also holds the following identity with the incompatibility operator \n",
    "\\begin{align*}\n",
    "\\mathrm{div}_g\\mathrm{div}_g(\\mathbb{S}_{g}\\sigma)\\omega(g) = -\\mathrm{inc}_g(\\sigma)\\omega(g),\\qquad \\mathrm{inc}_g:=\\mathrm{curl}_g\\mathrm{curl}_g.\n",
    "\\end{align*}\n",
    "In the Euclidean case, it would become $\\mathrm{div}\\mathrm{div}(\\mathbb{S}\\sigma)$. Using $\\sigma_h\\in \\mathrm{Reg}_h^k$ in the Regge space, $\\mathbb{S}\\sigma_h$ becomes normal-normal continuous fitting to the Hellan-Herrmann-Johnson method [Comodi. *The Hellan-Herrmann-Johnson method: some new error estimates and postprocessing*. *Math. Comp.* (1989)](https://doi.org/10.1090/S0025-5718-1989-0946601-7), which entails a well-defined distributional version reading for $u_h\\in\\mathrm{Lag}^{k+1}$ Lagrange finite elements as test functions\n",
    "\\begin{align*}\n",
    "\\langle \\mathrm{div}\\mathrm{div}(\\mathbb{S}\\sigma_h),u_h\\rangle &= \\sum_{T\\in\\mathcal{T}}\\int_T \\mathbb{S}\\sigma_h : \\nabla^2 u_h - \\int_{\\partial T} (\\mathbb{S}\\sigma_h)_{nn}\\partial_nu_h = \\sum_{T\\in\\mathcal{T}}\\int_T \\mathbb{S}\\sigma_h : \\nabla^2 u_h + \\int_{\\partial T} \\sigma_h(t,t)\\partial_nu_h\\\\\n",
    "&= \\sum_{T\\in\\mathcal{T}}\\int_T \\mathrm{div}\\mathrm{div}(\\mathbb{S}\\sigma_h) u_h + \\int_{\\partial T}\\mathrm{div}(\\mathbb{S}\\sigma)_n\\,u_h- (\\mathbb{S}\\sigma_h)_{nt}\\partial_tu_h\\\\\n",
    "&= \\sum_{T\\in\\mathcal{T}}\\int_T \\mathrm{div}\\mathrm{div}(\\mathbb{S}\\sigma_h) u_h + \\int_{\\partial T} u_h\\,\\left(\\mathrm{div}(\\mathbb{S}\\sigma)_n+\\partial_t(\\mathbb{S}\\sigma_h(n,t))\\right) + \\sum_{V}[\\![\\sigma_h(n,t)]\\!]_Vu_h(V),\n",
    "\\end{align*}\n",
    "where $[\\![\\sigma_h(n,t)]\\!]_V = \\sum_{T\\supset V}[\\![\\sigma_h(n,t)]\\!]^T_V$ and $[\\![\\sigma_h(n,t)]\\!]^T_V$ measures the jump of $\\sigma_h(n,t)$ evaluated at the triangle $T$ between the two edges of the triangle attached to the vertex $V$, i.e., $[\\![\\sigma_h(n,t)]\\!]^T_V = \\sigma_h(n,t)|_{E_1} -\\sigma_h(n,t)|_{E_2}$.\n",
    "\n",
    "Gawlik recognized that the HHJ method could be extended to the covariant case, which motivated to define for $g\\in \\mathrm{Reg}_h^k$ the following covariant version\n",
    "\\begin{align*}\n",
    "\\langle \\mathrm{div}_g\\mathrm{div}_g(\\mathbb{S}_g\\sigma_h)\\omega(g),u_h\\rangle &= \\sum_{T\\in\\mathcal{T}}\\int_T g(\\mathbb{S}_g\\sigma_h, \\nabla_g^2 u_h)\\,\\omega|_T(g) - \\int_{\\partial T} (\\mathbb{S}_g\\sigma_h)_{nn}\\partial_nu_h\\,\\omega_{\\partial T}(g)=\\sum_{T\\in\\mathcal{T}}\\int_T g(\\mathbb{S}_g\\sigma_h, \\nabla_g^2 u_h)\\,\\omega|_T(g) + \\int_{\\partial T} \\sigma_h(t,t)\\partial_nu_h\\,\\omega_{\\partial T}(g)\\\\\n",
    "&= \\sum_{T\\in\\mathcal{T}}\\int_T \\mathrm{div}\\mathrm{div}(\\mathbb{S}\\sigma_h) u_h + \\int_{\\partial T} \\mathrm{div}_g(\\mathbb{S}_g\\sigma)(n)\\,u_h-(\\mathbb{S}\\sigma_h)_{nt}\\partial_tu_h\\\\\n",
    "&= \\sum_{T\\in\\mathcal{T}}\\int_T \\mathrm{div}\\mathrm{div}(\\mathbb{S}\\sigma_h) u_h + \\int_{\\partial T}u_h\\, \\left(\\mathrm{div}_g(\\mathbb{S}_g\\sigma)(n)+\\partial_t(\\mathbb{S}\\sigma_h(n,t))\\right) + \\sum_{V}[\\![\\sigma_h(n,t)]\\!]_Vu_h(V).\n",
    "\\end{align*}\n",
    "It turns out that taking the variations of the jump of the geodesic curvature and the angle defect [Berchenko-Kogan, Gawlik. *Finite Element Approximation of the Levi-Civita Connection and Its Curvature in Two Dimensions*, *Found Comput Math* (2022)](https://doi.org/10.1007/s10208-022-09597-1) exactly match the terms of the covariant HHJ method\n",
    "\\begin{align*}\n",
    "D_g(\\kappa_g\\,\\omega_E)[\\sigma]=\\frac{1}{2}\\left( (\\mathrm{div}_g(\\mathbb{S}_g(\\sigma)))(n)+\\partial_{t}(\\sigma(n,t)) \\right)\\,\\omega_E,\\qquad D_g(\\sphericalangle_V)[\\sigma]=\\frac{1}{2}\\sum_{T\\supset V}[\\![\\sigma(n,t)]\\!]_V^T\n",
    "\\end{align*}\n",
    "proving that\n",
    "\\begin{align*}\n",
    "D_g \\widetilde{K\\omega}(u_h)[\\sigma]=\\frac{1}{2}b_h(g;\\sigma_h,u_h),\\qquad b_h(g;\\sigma_h,u_h) =\\sum_{T\\in\\mathcal{T}}\\int_T g(\\mathbb{S}_g\\sigma_h, \\nabla_g^2 u_h)\\,\\omega|_T(g) - \\int_{\\partial T} \\sigma_h(t,t)\\partial_nu_h\\,\\omega_{\\partial T}(g)\n",
    "\\end{align*}\n",
    "Noting that the Gauss curvature for the Euclidean metric is zero, we obtain with the main theorem of integration and differentiation the following integral representation of the Gauss curvature\n",
    "\\begin{align*}\n",
    "\\widetilde{K\\omega}(u_h)= \\frac{1}{2}\\int_0^1b_h(g(t);\\sigma,u_h)\\,dt,\\qquad g(t) = \\delta+t(g_h-\\delta),\\quad \\sigma=\\dot{g}(t) = g_h-\\delta,\n",
    "\\end{align*}\n",
    "where $\\delta$ denotes the Euclidean metric.\n",
    "By consistency, the exact Gauss curvature from the smooth, exact metric also fulfils the integral representation. Thus, the error has the following integral representation\n",
    "\\begin{align*}\n",
    "\\langle \\widetilde{K\\omega}(g_h)-K\\omega,u_h\\rangle = \\frac{1}{2}\\int_0^1b_h(g(t);\\sigma,u_h)\\,dt,\\qquad g(t) = g+t(g_h-g),\\quad \\sigma=\\dot{g}(t) = g_h-g\n",
    "\\end{align*}\n",
    "The integrand is linear enough to perform numerical analysis to prove the observed convergence rates. For the details of the (technical) proofs, we refer to the above-mentioned literature."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "606eac78",
   "metadata": {},
   "source": [
    "## Improved convergence rates by the canonical Regge interpolant\n",
    "When using the canonical Regge interpolation operator to approximate the metric tensor, i.e. $g_h=\\mathcal{I}^{\\mathrm{Reg}^k}g$, one can achieve one or even two ordesr of convergence more. For $g_h=\\mathcal{I}^{\\mathrm{Reg}^{k-1}}g$ and $K_h\\in\\mathrm{Lag}^k$ there holds the following convergence rates\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\|K_h-K\\|_{L^2}\\le C h^{k},\\qquad \\|K_h-K\\|_{H^{-1}}\\le C h^{k+1}.\n",
    "\\end{align*}\n",
    "$$\n",
    "Using stronger norms is also possible and reduces the rates as usual. The proof's strategy relies on exploiting the definition of the canonical Regge interpolant that the first $k+1$ tangential-tangential moments at the edges (starting with the zeroth-order moment) and, for higher-order, some internal moments of $g$ and $g_h$ coincide. This allows the extraction of one additional order of convergence. The pretty technical proof can be found in [Gopalakrishnan, Neunteufel, Schöberl, Wardetzky. *Analysis of curvature approximations via covariant curl and incompatibility for Regge metrics*, *The SMAI Journal of computational mathematics* (2023).](https://doi.org/10.5802/smai-jcm.98).\n",
    "\n",
    "\n",
    "\n",
    "By reducing the polynomial order of $K_h$ to coincide with the metric's order, i.e. $K_h\\in \\mathrm{Lag}^k$ and $g_h=\\mathcal{I}^{\\mathrm{Reg}^k}g$, under the assumption that $k>1$, yields yet another improvement of the convergence rate\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+2},\\qquad k>1.\n",
    "\\end{align*}\n",
    "$$\n",
    "In fact, one obtains the optimal convergence possible in the $L^2$-norm. Again, the key ingredient is the definition of the canonical Regge interpolant. The details can be found in [Gopalakrishnan, Neunteufel, Schöberl, Wardetzky. *On the improved convergence of lifted distributional Gauss curvature from Regge elements*, *Results in Applied Mathematics* (2024).](https://doi.org/10.1016/j.rinam.2024.100511). Note that the convergence rates highly depend on the chosen polynomial approximation spaces for $g_h$ and $K_h$, and whether the canonical Regge interpolation operator is used or something else."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43e1528e",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}