{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3561e8aa",
   "metadata": {},
   "source": [
    "# Distributional Einstein curvature\n",
    "This notebook discusses the definition and convergence of the distributional Einstein tensor in dimensions three and higher (in one and two dimensions, the Einstein tensor is zero). For a detailed derivation and analysis, we refer to [Gawlik, Neunteufel. *Finite element approximation of the Einstein tensor*, *IMA Journal of Numerical Analysis* (2025).](https://doi.org/10.1093/imanum/draf004).\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "## Einstein tensor as a nonlinear distribution\n",
    "The Einstein tensor reads with the Ricci curvature tensor $\\mathrm{Ric}$ and scalar curvature $S$\n",
    "\\begin{align*}\n",
    "G = \\mathrm{Ric}-\\frac{1}{2}S\\,g.\n",
    "\\end{align*}\n",
    "One important application of the Einstein tensor for pseudo-Riemannian manifolds is general relativity. The Einstein field equations $G=8\\pi\\mathfrak{T}$, where $\\mathfrak{T}$ is the energy-stress tensor, describe cosmologic phenomena such as black holes.\n",
    "\n",
    "\n",
    "Motivated by the integral representation of the [distributional scalar curvature](distributional_scalar_curvature.ipynb), we define the distributional (densitized) Einstein tensor acting on a sufficiently smooth matrix-valued test function $\\varphi$ as\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\widetilde{G\\omega}(\\varphi) :=\\sum_{T\\in\\mathcal{T}}\\int_T\\langle G|_T,\\varphi\\rangle\\,\\omega_T+\\sum_{F\\in\\mathcal{F}}\\int_F\\langle[\\![\\bar{I\\!I}]\\!]_F,\\varphi|_F\\rangle\\,\\omega_F-\\sum_{E\\in\\mathcal{E}}\\int_E\\langle\\Theta_Eg|_E,\\varphi|_E\\rangle\\,\\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$. $\\bar{I\\!I}=I\\!I-H\\,g$ is the trace-reversed second fundamental form with mean curvature $H$\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: Different sign conventions exist again). $\\Theta_E$ is the angle defect. As we have seen in the previous notebook, the concept of the angle defect directly extends from 2D to any dimension. \n",
    "\n",
    "The jump of the trace-reversed second fundamental form at facets is well-known in physics literature as the Israel junction conditions [Israel. *Singular hypersurfaces and thin shells in general relativity*, *Nuovo Cimento B (1965-1970)* 1966.](https://doi.org/10.1007/BF02710419) verifying our distributional facet term."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91fc6eab",
   "metadata": {},
   "source": [
    "## Integral representation and convergence\n",
    "The distributional EInstein tensor entails an integral representation, which is more involved than the scalar curvature.\n",
    "Taking the variation of the densitized distributional Einstein tensor in the metric in direction $\\sigma$ yields\n",
    "$$\n",
    "\\begin{align*}\n",
    "D_g\\widetilde{G\\omega}(\\rho)[\\sigma] = B_h(g;\\sigma,\\rho)-A_h(g;\\sigma,\\rho)\n",
    "\\end{align*}\n",
    "$$\n",
    "with\n",
    "$$\n",
    "\\begin{align*}\n",
    "2B_h(g;\\sigma,\\rho)&= \\sum_T\\int_T\\langle 2\\mathrm{ein}\\sigma,\\rho\\rangle\\,\\omega_T+\\sum_{F}\\int_F\\langle[\\![\\mathbb{S}_F(\\sigma(n,n)I\\!I+\\nabla_n\\sigma-(\\nabla^F\\sigma)(n,\\cdot)-\\nabla^F(\\sigma(n,\\cdot)))]\\!],\\rho|_F\\rangle\\,\\omega_F\\\\\n",
    "&\\qquad-\\sum_E\\int_E\\sum_{F\\supset E}[\\![\\sigma(n,\\nu)]\\!]_F\\mathrm{tr}(\\rho|_E)\\,\\omega_E,\\\\\n",
    "2A_h(g;\\sigma,\\rho)&= \\sum_T\\int_T\\left(2\\sigma:\\mathfrak{R}:\\rho +\\langle\\mathrm{Ric},\\sigma \\rangle\\,\\mathrm{tr}\\rho+\\langle\\mathrm{Ric},\\rho \\rangle\\,\\mathrm{tr}\\sigma + S\\langle J\\sigma,\\rho\\rangle-2\\sigma:\\mathrm{Ric}:\\rho\\right)\\,\\omega_T\\\\\n",
    "&\\qquad+\\sum_F\\int_F\\left(-3(\\sigma|_F):[\\![\\bar{I\\!I}]\\!]:(\\rho|_F)+\\langle [\\![\\bar{I\\!I}]\\!],\\sigma|_F\\rangle\\,\\mathrm{tr}(\\rho_F)+\\langle [\\![\\bar{I\\!I}]\\!],\\rho|_F\\rangle\\,\\mathrm{tr}(\\sigma|_F)-[\\![H]\\!]\\langle\\mathbb{S}_F\\sigma,\\rho|_F\\rangle\\right)\\omega_F\\\\\n",
    "&\\qquad+\\sum_{E}\\int_E\\left(\\langle\\Theta_E\\sigma|_E,\\rho|_E\\rangle-\\Theta_E\\mathrm{tr}(\\sigma|_E)\\,\\mathrm{tr}(\\rho|_E)\\right)\\,\\omega_S.\n",
    "\\end{align*}\n",
    "$$\n",
    "Here, all differential operators are covariant ones with respect to the metric $g$, and $\\nu$ denotes the co-normal vector being perpendicular to the co-dimension 2 boundary $E$ pointing from $E$ into the adjacent facets $F$ of the element $T$. $\\mathrm{ein}$ is the (covariant) linearized Einstein operator defined by\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\mathrm{ein} = J\\,\\mathrm{def}\\,\\mathrm{div}\\,J-\\frac{1}{2}\\Delta\\,J,\\qquad \\mathrm{def} = \\mathrm{sym}\\,\\nabla,\\quad J\\sigma =\\sigma -\\frac{1}{2}\\mathrm{tr}\\sigma\\,g\n",
    "\\end{align*}\n",
    "$$\n",
    "and $\\mathbb{S}_F\\sigma=\\sigma|_F-\\mathrm{tr}(\\sigma|_F)g|_F$. Note that $A(g;\\sigma,\\rho)$ and $B(g;\\sigma,\\rho)$ become identically zero in two dimensions. $B_h(g;\\sigma,\\rho)$ can be interpreted as a distributional (covariant) extension of the linearized Einstein operator $\\mathrm{ein}$ for tangential-tangential continuous matrices. In three dimensions, the linearized Einstein operator coincides with the incompatibility operator by up to a factor of $2$.\n",
    "\n",
    "The convergence rates in the $H^{-2}$ norm for $k\\ge 0$ and $g_h\\in \\mathrm{Reg}^k$ coincide with those of the distributional scalar curvature\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\|\\widetilde{G\\omega}(g_h)-(G\\omega)(g)\\|_{H^{-2}}\\le C\n",
    "h^{\\min\\{k+1,2k\\}}.\n",
    "\\end{align*}\n",
    "$$\n",
    "Again, we cannot expect convergence in the lowest-order case of a piecewise constant Regge metric. 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."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43e1528e",
   "metadata": {},
   "source": [
    "## Numerical example\n",
    "In the following example, we compute the distributional Einstein tensor $\\widetilde{G\\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=12, ny=12, nz=12, mapping=mapping)\n",
    "\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",
    "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",
    "def q(X):\n",
    "    return X**2 * (X**2 - 3) ** 2\n",
    "\n",
    "\n",
    "# exact Ricci curvature tensor\n",
    "Ricci11 = (\n",
    "    -9\n",
    "    * (x**2 - 1)\n",
    "    * ((y**2 + z**2 - 2) * (q(x) + 9) + (z**2 - 1) * q(y) + q(z) * (y**2 - 1))\n",
    "    / (9 + q(x) + q(y) + q(z)) ** 2\n",
    ")\n",
    "Ricci22 = (\n",
    "    -9\n",
    "    * (y**2 - 1)\n",
    "    * ((x**2 + z**2 - 2) * (q(y) + 9) + (z**2 - 1) * q(x) + q(z) * (x**2 - 1))\n",
    "    / (9 + q(x) + q(y) + q(z)) ** 2\n",
    ")\n",
    "Ricci33 = (\n",
    "    -9\n",
    "    * (z**2 - 1)\n",
    "    * ((y**2 + x**2 - 2) * (q(z) + 9) + (x**2 - 1) * q(y) + q(x) * (y**2 - 1))\n",
    "    / (9 + q(x) + q(y) + q(z)) ** 2\n",
    ")\n",
    "Ricci12 = (\n",
    "    -9\n",
    "    * (y**2 - 3)\n",
    "    * y\n",
    "    * (x**2 - 3)\n",
    "    * x\n",
    "    * (x**2 - 1)\n",
    "    * (y**2 - 1)\n",
    "    / (9 + q(x) + q(y) + q(z)) ** 2\n",
    ")\n",
    "Ricci13 = (\n",
    "    -9\n",
    "    * (z**2 - 3)\n",
    "    * z\n",
    "    * (x**2 - 3)\n",
    "    * x\n",
    "    * (x**2 - 1)\n",
    "    * (z**2 - 1)\n",
    "    / (9 + q(x) + q(y) + q(z)) ** 2\n",
    ")\n",
    "Ricci23 = (\n",
    "    -9\n",
    "    * (y**2 - 3)\n",
    "    * y\n",
    "    * (z**2 - 3)\n",
    "    * z\n",
    "    * (z**2 - 1)\n",
    "    * (y**2 - 1)\n",
    "    / (9 + q(x) + q(y) + q(z)) ** 2\n",
    ")\n",
    "\n",
    "Ricci_ex = -CF(\n",
    "    (Ricci11, Ricci12, Ricci13, Ricci12, Ricci22, Ricci23, Ricci13, Ricci23, Ricci33),\n",
    "    dims=(3, 3),\n",
    ")\n",
    "\n",
    "# exact scalar curvature\n",
    "S_ex = (\n",
    "    -(\n",
    "        2 * ((1 - z**2) * y**2 + z**2 - 1) * x**12\n",
    "        + 24 * ((-1 + z**2) * y**2 - z**2 + 1) * x**10\n",
    "        + (\n",
    "            2 * (1 - z**2) * y**6\n",
    "            + 12 * (-1 + z**2) * y**4\n",
    "            + (108 - 2 * z**6 - 144 * z**2 + 12 * z**4) * y**2\n",
    "            - 72\n",
    "            - 12 * z**4\n",
    "            + 108 * z**2\n",
    "            + 2 * z**6\n",
    "        )\n",
    "        * x**8\n",
    "        + (\n",
    "            2 * (1 - z**2) * y**8\n",
    "            + 28 * (-1 + z**2) * y**6\n",
    "            + 114 * (1 - z**2) * y**4\n",
    "            + (468 * z**2 + 28 * z**6 - 2 * z**8 - 198 - 114 * z**4) * y**2\n",
    "            - 72\n",
    "            + 114 * z**4\n",
    "            + 2 * z**8\n",
    "            - 28 * z**6\n",
    "            - 198 * z**2\n",
    "        )\n",
    "        * x**6\n",
    "        + (\n",
    "            (-12 + 12 * z**2) * y**8\n",
    "            + (114 - 114 * z**2) * y**6\n",
    "            + (-360 + 360 * z**2) * y**4\n",
    "            + (360 * z**4 - 114 * z**6 + 12 * z**8 + 54 - 702 * z**2) * y**2\n",
    "            + 594\n",
    "            + 114 * z**6\n",
    "            - 12 * z**8\n",
    "            - 360 * z**4\n",
    "            + 54 * z**2\n",
    "        )\n",
    "        * x**4\n",
    "        + (\n",
    "            (2 - 2 * z**2) * y**12\n",
    "            + (-24 + 24 * z**2) * y**10\n",
    "            + (108 - 2 * z**6 - 144 * z**2 + 12 * z**4) * y**8\n",
    "            + (468 * z**2 + 28 * z**6 - 2 * z**8 - 198 - 114 * z**4) * y**6\n",
    "            + (360 * z**4 - 114 * z**6 + 12 * z**8 + 54 - 702 * z**2) * y**4\n",
    "            + (486 + 468 * z**6 - 2 * z**12 + 24 * z**10 - 144 * z**8 - 702 * z**4)\n",
    "            * y**2\n",
    "            + 108 * z**8\n",
    "            - 324\n",
    "            - 24 * z**10\n",
    "            + 2 * z**12\n",
    "            - 198 * z**6\n",
    "            + 54 * z**4\n",
    "            + 486 * z**2\n",
    "        )\n",
    "        * x**2\n",
    "        + (2 * z**2 - 2) * y**12\n",
    "        + (-24 * z**2 + 24) * y**10\n",
    "        + (-72 - 12 * z**4 + 108 * z**2 + 2 * z**6) * y**8\n",
    "        + (-72 + 114 * z**4 + 2 * z**8 - 28 * z**6 - 198 * z**2) * y**6\n",
    "        + (594 + 114 * z**6 - 12 * z**8 - 360 * z**4 + 54 * z**2) * y**4\n",
    "        + (\n",
    "            108 * z**8\n",
    "            - 324\n",
    "            - 24 * z**10\n",
    "            + 2 * z**12\n",
    "            - 198 * z**6\n",
    "            + 54 * z**4\n",
    "            + 486 * z**2\n",
    "        )\n",
    "        * y**2\n",
    "        - 486\n",
    "        - 72 * z**8\n",
    "        + 24 * z**10\n",
    "        - 2 * z**12\n",
    "        - 72 * z**6\n",
    "        + 594 * z**4\n",
    "        - 324 * z**2\n",
    "    )\n",
    "    / (\n",
    "        1\n",
    "        + 1 / 9 * x**6\n",
    "        + 1 / 9 * y**6\n",
    "        + 1 / 9 * z**6\n",
    "        - 2 / 3 * y**4\n",
    "        + y**2\n",
    "        - 2 / 3 * x**4\n",
    "        + x**2\n",
    "        - 2 / 3 * z**4\n",
    "        + z**2\n",
    "    )\n",
    "    / (\n",
    "        x**6\n",
    "        + y**6\n",
    "        + z**6\n",
    "        - 6 * x**4\n",
    "        - 6 * y**4\n",
    "        - 6 * z**4\n",
    "        + 9 * x**2\n",
    "        + 9 * y**2\n",
    "        + 9 * z**2\n",
    "        + 9\n",
    "    )\n",
    "    ** 2\n",
    ")\n",
    "\n",
    "# exact Einstein tensor\n",
    "Einstein_ex = Ricci_ex - 1 / 2 * S_ex * G_ex\n",
    "clipping = {\"pnt\": (0, 0, 0.01), \"z\": -1, \"function\": True}\n",
    "Draw(Norm(Einstein_ex), mesh, clipping=clipping);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d35a8ac3",
   "metadata": {},
   "source": [
    "We compute the $H^{-2}$-norm of the error $\\|\\widetilde{G\\omega}-G\\omega\\|_{H^{-2}}$ by solving for each component the following fourth-order problem $\\mathrm{div}(\\mathrm{div}(\\nabla^2 u_{ij}))= (\\widetilde{G\\omega}-G\\omega)_{ij}$. Then $\\|u\\|_{H^2}\\approx \\|\\widetilde{G\\omega}-G\\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_{ij}$ 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": [
    "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",
    "    apre = Preconditioner(a, \"bddc\")\n",
    "    a.Assemble()\n",
    "\n",
    "    invS = CGSolver(a.mat, apre, printrates=\"\\r\", maxiter=600)\n",
    "    if a.condense:\n",
    "        ext = IdentityMatrix() + a.harmonic_extension\n",
    "        inv = a.inner_solve + ext @ invS @ ext.T\n",
    "    else:\n",
    "        inv = invS\n",
    "\n",
    "    # use MultiVector to solve for all right-hand sides at once\n",
    "    rhs = MultiVector(fes.ndof, 6, False)\n",
    "    gf_u_i = GridFunction(fes**6)\n",
    "\n",
    "    for i in range(6):\n",
    "        j = i if i < 3 else (i + 1 if i < 5 else 8)\n",
    "        f = LinearForm(\n",
    "            func_vol[j] * v * dx\n",
    "            + func_bnd[j] * v * dx(element_boundary=True)\n",
    "            + func_bbnd[j] * v * dx(element_vb=BBND)\n",
    "        )\n",
    "        f.Assemble()\n",
    "        rhs[i].data = f.vec\n",
    "\n",
    "    gfu_multivec = (inv * rhs).Evaluate()\n",
    "\n",
    "    for i in range(6):\n",
    "        gf_u_i.components[i].vec.data = gfu_multivec[i]\n",
    "\n",
    "    index_symmetric_matrix = [0, 1, 2, 1, 3, 4, 2, 4, 5]\n",
    "    # u\n",
    "    gfu = CF(\n",
    "        tuple([gf_u_i.components[i].components[0] for i in index_symmetric_matrix]),\n",
    "        dims=(3, 3),\n",
    "    )\n",
    "    # Grad(u)\n",
    "    grad_gfu = CF(\n",
    "        tuple(\n",
    "            [Grad(gf_u_i.components[i].components[0]) for i in index_symmetric_matrix]\n",
    "        ),\n",
    "        dims=(3, 3, 3),\n",
    "    )\n",
    "    # sigma = hesse u\n",
    "    gf_sigma = CF(\n",
    "        tuple([gf_u_i.components[i].components[1]] for i in index_symmetric_matrix),\n",
    "        dims=(3, 3, 3, 3),\n",
    "    )\n",
    "    err = sqrt(\n",
    "        Integrate(\n",
    "            InnerProduct(gfu, gfu)\n",
    "            + InnerProduct(grad_gfu, grad_gfu)\n",
    "            + InnerProduct(gf_sigma, gf_sigma),\n",
    "            mesh,\n",
    "        )\n",
    "    )\n",
    "    return err"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74acccd7",
   "metadata": {},
   "source": [
    "The following function computes the $H^2$-error of the distributional Einstein tensor for a given metric and polynomial order. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5f0e469",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ComputeEinsteinCurvature(mesh, G, order):\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",
    "    n = specialcf.normal(mesh.dim)\n",
    "\n",
    "    # Regge metric\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",
    "    # co-dimension 2 term\n",
    "    n1g = 1 / sqrt(Inv(gfG) * n1 * n1) * Inv(gfG) * n1\n",
    "    n2g = 1 / sqrt(Inv(gfG) * n2 * n2) * Inv(gfG) * n2\n",
    "    tEg = 1 / sqrt(gfG * tE * tE) * tE\n",
    "    func_bbnd = (\n",
    "        -vol_bbnd * (acos(n1 * n2) - acos(gfG * n1g * n2g)) * OuterProduct(tEg, tEg)\n",
    "    )\n",
    "\n",
    "    # co-dimension 1 facet term\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 = vol_bnd * P * (SFF - H * gfG) * P\n",
    "\n",
    "    # volume term\n",
    "    func_vol = vol * Inv(gfG) * gfG.Operator(\"Einstein\") * Inv(gfG) - sqrt(\n",
    "        Det(G_ex)\n",
    "    ) * (Inv(G_ex) * Einstein_ex * Inv(G_ex))\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 Einstein tensor 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(4):\n",
    "        mesh = GenerateRandomMesh(k=i, par=3)\n",
    "        errm2, dof = ComputeEinsteinCurvature(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{G\\omega}-G_{\\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": "fb563e6a",
   "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
}