{ "cells": [ { "cell_type": "markdown", "id": "192c960a", "metadata": {}, "source": [ "# Linearization of curvature quantities\n", "\n", "In the numerical analysis of distributional curvatures, their linearization around a given metric is essential, i.e., how the curvature quantity transforms when the underlying metric changes. These linearizations involve covariant differential operators, which we are going to investigate in this notebook by testing the formulas numerically by means of a Taylor test. To this end we compute the quantity $F=F(g)$ and a time-dependent $g=g(t)$ and denote $\\sigma = \\dot{g}$\n", "$$\n", "F(g+t\\sigma)-F(g)-t\\, \\dot{F}(g),\n", "$$\n", "where $\\dot{F}(g):= \\frac{d}{dt} F(g):=(DF)(g)[\\sigma]$ denotes the directional derivative in direction $\\sigma$ and test if the expression converges quadratically in $t$ for $t\\to 0$." ] }, { "cell_type": "markdown", "id": "66eadd48", "metadata": {}, "source": [ "## 2D" ] }, { "cell_type": "code", "execution_count": null, "id": "1f7ff184", "metadata": {}, "outputs": [], "source": [ "from math import log2\n", "from ngsolve import *\n", "from netgen.occ import *\n", "import ngsdiffgeo as dg\n", "from ngsolve.meshes import MakeStructured2DMesh\n", "\n", "\n", "mesh = MakeStructured2DMesh(False, nx=6, ny=6)\n", "peak = 0.5 * x**2 - 1 / 12 * x**4 + 0.5 * y**2 - 1 / 12 * y**4\n", "F = CF((1, 0, 0, 1, peak.Diff(x), peak.Diff(y)), dims=(3, 2))\n", "Gex = F.trans * F\n", "\n", "cfsigma = dg.TensorField(\n", " Sym(\n", " CF(\n", " (\n", " 1 + 10 * x * y**3 - x**2,\n", " 0.2 * y**4 * x - y,\n", " 0.2 * sin(x * y),\n", " 1 + cos(x) * y**2,\n", " ),\n", " dims=(2, 2),\n", " )\n", " ),\n", " \"11\",\n", ")\n", "\n", "order = 3\n", "\n", "gfG = GridFunction(HCurlCurl(mesh, order=order))\n", "gfsigma = GridFunction(HCurlCurl(mesh, order=order))\n", "\n", "with TaskManager():\n", " gfG.Set(Gex)\n", " gfsigma.Set(cfsigma)\n", "\n", "sigma = dg.TensorField(gfsigma, \"11\")\n", "\n", "mf = dg.RiemannianManifold(gfG)" ] }, { "cell_type": "markdown", "id": "fcff9851", "metadata": {}, "source": [ "We start by verifying the variations of the volume forms\n", "\n", "$$\n", "\\frac{d}{dt}(\\omega_T) = 0.5\\,\\mathrm{tr}(\\sigma)\\,\\omega_T,\\qquad \\frac{d}{dt}(\\omega_F) = 0.5\\,\\mathrm{tr}(\\sigma|_F)\\,\\omega_F\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "id": "f55ac670", "metadata": {}, "outputs": [], "source": [ "t = 1 / 8\n", "gfGt = GridFunction(gfsigma.space)\n", "errold = None\n", "err = None\n", "\n", "omega_T = mf.VolumeForm(VOL)\n", "omega_F = mf.VolumeForm(BND)\n", "\n", "print(\"omega_T:\")\n", "with TaskManager():\n", " term_0 = omega_T\n", "\n", " variation = 0.5 * mf.Trace(sigma) * omega_T\n", "\n", " for m in range(8):\n", " t /= 2\n", " gfGt.vec.data = gfG.vec + t * gfsigma.vec\n", " mf_t = dg.RiemannianManifold(gfGt)\n", " term_t = mf_t.VolumeForm(VOL)\n", "\n", " errold = err\n", " err = sqrt(Integrate((term_t - term_0 - t * variation) ** 2, mesh))\n", " print(\n", " f\"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}\"\n", " )\n", "\n", "\n", "print(\"\\nomega_F:\")\n", "t = 1 / 8\n", "errold = None\n", "err = None\n", "\n", "with TaskManager():\n", " term_0 = omega_F\n", "\n", " variation = 0.5 * mf.Trace(sigma, BND) * omega_F\n", "\n", " for m in range(8):\n", " t /= 2\n", " gfGt.vec.data = gfG.vec + t * gfsigma.vec\n", " mf_t = dg.RiemannianManifold(gfGt)\n", " term_t = mf_t.VolumeForm(BND)\n", "\n", " errold = err\n", "\n", " err = sqrt(\n", " Integrate(\n", " (term_t - term_0 - t * variation) ** 2 * dx(element_boundary=True), mesh\n", " )\n", " )\n", " print(\n", " f\"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}\"\n", " )" ] }, { "cell_type": "markdown", "id": "094fd30d", "metadata": {}, "source": [ "Next, we look at the Gauss curvature $K$ and the Gauss curvature multiplied by the volume form $K\\,\\omega_T$. The variations are\n", "\n", "$$\n", "\\frac{d}{dt}K = 0.5(\\mathrm{div}_g\\mathrm{div}_g(\\sigma)-\\Delta_g\\mathrm{tr}(\\sigma)-g(\\sigma,\\mathrm{Ric})),\\qquad \\frac{d}{dt}(K\\,\\omega_T) = 0.5\\,\\mathrm{div}_g\\mathrm{div}_g(\\mathbb{S}_g\\sigma)\\,\\omega_T,\n", "$$\n", "where $\\mathrm{div}_g$ is the covariant divergence of vector/tensor-fields, $\\Delta_g$ the Laplace-Beltrami operator, and $\\mathbb{S}_g\\sigma:=\\sigma-\\mathrm{tr}(\\sigma)\\,g$." ] }, { "cell_type": "code", "execution_count": null, "id": "989e4a95", "metadata": {}, "outputs": [], "source": [ "print(\"Gauss curvature:\")\n", "t = 1 / 8\n", "gfGt = GridFunction(gfsigma.space)\n", "errold = None\n", "err = None\n", "\n", "with TaskManager():\n", " term_0 = mf.Gauss\n", "\n", " variation = 0.5 * (\n", " mf.CovDiv(mf.CovDiv(sigma))\n", " - mf.CovLaplace(mf.Trace(sigma))\n", " - mf.InnerProduct(sigma, mf.Ricci)\n", " )\n", "\n", " for m in range(8):\n", " t /= 2\n", " gfGt.vec.data = gfG.vec + t * gfsigma.vec\n", " mf_t = dg.RiemannianManifold(gfGt)\n", " term_t = mf_t.Gauss\n", "\n", " errold = err\n", " err = sqrt(Integrate((term_t - term_0 - t * variation) ** 2, mesh))\n", " print(\n", " f\"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}\"\n", " )\n", "\n", "\n", "print(\"\\nK omega_T:\")\n", "t = 1 / 8\n", "gfGt = GridFunction(gfsigma.space)\n", "errold = None\n", "err = None\n", "\n", "with TaskManager():\n", " term_0 = mf.Gauss * omega_T\n", "\n", " variation = 0.5 * mf.CovDiv(mf.CovDiv(mf.S(sigma))) * omega_T\n", "\n", " for m in range(8):\n", " t /= 2\n", " gfGt.vec.data = gfG.vec + t * gfsigma.vec\n", " mf_t = dg.RiemannianManifold(gfGt)\n", " term_t = mf_t.Gauss * mf_t.VolumeForm(VOL)\n", "\n", " errold = err\n", " err = sqrt(Integrate((term_t - term_0 - t * variation) ** 2, mesh))\n", " print(\n", " f\"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}\"\n", " )" ] }, { "cell_type": "markdown", "id": "eca74a42", "metadata": {}, "source": [ "For the geodesic curvature $\\kappa$ at edges we have that\n", "$$\n", "\\frac{d}{dt}(\\kappa\\,\\omega_F) = 0.5\\Big(g(\\mathrm{div}_g(\\mathbb{S}_g\\sigma), n_g) + g\\big(\\nabla_g(\\sigma(n_g,t_g)),t_g\\big) \\Big)\\,\\omega_F,\n", "$$\n", "where $n_g$ and $t_g$ are the $g$-normal and tangential vectors." ] }, { "cell_type": "code", "execution_count": null, "id": "506267f4", "metadata": {}, "outputs": [], "source": [ "print(\"kappa omega_F:\")\n", "t = 1 / 8\n", "gfGt = GridFunction(gfsigma.space)\n", "errold = None\n", "err = None\n", "\n", "tv = specialcf.tangential(mesh.dim)\n", "nv = specialcf.normal(mesh.dim)\n", "\n", "nv_g = dg.VectorField(mf.normal)\n", "tv_g = mf.tangent\n", "\n", "with TaskManager():\n", " term_0 = mf.GeodesicCurvature * omega_F\n", "\n", " variation = (\n", " 0.5\n", " * (\n", " mf.InnerProduct(mf.CovDiv(mf.S(sigma)), nv_g)\n", " + mf.InnerProduct(\n", " mf.CovDeriv(\n", " dg.ScalarField(mf.InnerProduct(sigma, dg.TensorProduct(nv_g, tv_g)))\n", " ),\n", " tv_g,\n", " )\n", " )\n", " * omega_F\n", " )\n", "\n", " for m in range(8):\n", " t /= 2\n", " gfGt.vec.data = gfG.vec + t * gfsigma.vec\n", " mf_t = dg.RiemannianManifold(gfGt)\n", " term_t = mf_t.GeodesicCurvature * mf_t.VolumeForm(BND)\n", "\n", " errold = err\n", " err = sqrt(\n", " Integrate(\n", " (term_t - term_0 - t * variation) ** 2 * dx(element_boundary=True), mesh\n", " )\n", " )\n", " print(\n", " f\"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}\"\n", " )" ] }, { "cell_type": "markdown", "id": "5d2ac3d8", "metadata": {}, "source": [ "## 3D" ] }, { "cell_type": "code", "execution_count": null, "id": "2afca643", "metadata": {}, "outputs": [], "source": [ "from math import log2\n", "from ngsolve import *\n", "from ngsolve.fem import Einsum\n", "from netgen.occ import *\n", "import ngsdiffgeo as dg\n", "from ngsolve.meshes import MakeStructured3DMesh\n", "\n", "mesh = MakeStructured3DMesh(False, nx=2, ny=2, nz=2)\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", "F = CF(\n", " (1, 0, 0, 0, 1, 0, 0, 0, 1, peak.Diff(x), peak.Diff(y), peak.Diff(z)),\n", " dims=(4, 3),\n", ")\n", "Gex = F.trans * F\n", "cfsigma = dg.TensorField(\n", " Sym(\n", " CF(\n", " (\n", " 1 + 10 * x * y**3 - x**2,\n", " 0.2 * y**4 * x - y,\n", " 0.2 * z**4 * x - z,\n", " 0.2 * sin(x * y),\n", " 1 + cos(x) * y**2,\n", " 0.2 * sin(x * z),\n", " 0.2 * z**4 * x - z,\n", " 0.2 * sin(x * z),\n", " 1 + cos(x) * z**2,\n", " ),\n", " dims=(3, 3),\n", " )\n", " ),\n", " \"11\",\n", ")\n", "order = 2\n", "\n", "gfG = GridFunction(HCurlCurl(mesh, order=order))\n", "gfsigma = GridFunction(HCurlCurl(mesh, order=order))\n", "\n", "with TaskManager():\n", " gfG.Set(Gex)\n", " gfsigma.Set(cfsigma)\n", "\n", "sigma = dg.TensorField(gfsigma, \"11\")\n", "\n", "mf = dg.RiemannianManifold(gfG)" ] }, { "cell_type": "markdown", "id": "489e2d11", "metadata": {}, "source": [ "In 3D, we also have the volume form on codimension 2 boundaries\n", "$$\n", "\\frac{d}{dt}(\\omega_E) = 0.5\\,\\mathrm{tr}(\\sigma|_E)\\,\\omega_E.\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "id": "6665c48c", "metadata": {}, "outputs": [], "source": [ "t = 1 / 8\n", "gfGt = GridFunction(gfsigma.space)\n", "errold = None\n", "err = None\n", "\n", "omega_T = mf.VolumeForm(VOL)\n", "omega_F = mf.VolumeForm(BND)\n", "omega_E = mf.VolumeForm(BBND)\n", "\n", "print(\"omega_E:\")\n", "t = 1 / 8\n", "errold = None\n", "err = None\n", "\n", "with TaskManager():\n", " term_0 = omega_E\n", "\n", " variation = 0.5 * mf.Trace(sigma, BBND) * omega_E\n", "\n", " for m in range(10):\n", " t /= 2\n", " gfGt.vec.data = gfG.vec + t * gfsigma.vec\n", " mf_t = dg.RiemannianManifold(gfGt)\n", " term_t = mf_t.VolumeForm(BBND)\n", "\n", " f = BilinearForm(gfsigma.space)\n", " f += Variation(\n", " InnerProduct(\n", " term_t - term_0 - t * variation, term_t - term_0 - t * variation\n", " )\n", " * dx(element_vb=BBND)\n", " )\n", " errold = err\n", " err = sqrt(f.Energy(gfsigma.vec))\n", " print(\n", " f\"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}\"\n", " )" ] }, { "cell_type": "markdown", "id": "b504a621", "metadata": {}, "source": [ "The scalar curvature coincides with the Gauss curvature in two dimensions up to a factor 2. In 3D, the linearization of the scalar curvature reads\n", "$$\n", "\\begin{align*}\n", "\\frac{d}{dt}S &= \\mathrm{div}_g\\mathrm{div}_g\\sigma-\\Delta_g\\mathrm{tr}\\sigma-g(\\mathrm{Ric},\\sigma),\\\\\n", "\\frac{d}{dt}(S\\,\\omega_T) &= \\big(\\mathrm{div}_g\\mathrm{div}_g(\\mathbb{S}_g\\sigma)-g(\\sigma,G)\\big)\\,\\omega_T,\n", "\\end{align*}\n", "$$\n", "where $\\Delta_g$ is the Laplace-Beltrami operator, $\\mathrm{Ric}$ the Ricci curvature tensor, and $G$ the Einstein tensor." ] }, { "cell_type": "code", "execution_count": null, "id": "81ae2901", "metadata": {}, "outputs": [], "source": [ "print(\"\\nscalar:\")\n", "t = 1 / 8\n", "gfGt = GridFunction(gfsigma.space)\n", "errold = None\n", "err = None\n", "\n", "with TaskManager():\n", " term_0 = mf.Scalar\n", " variation = (\n", " mf.CovDiv(mf.CovDiv(sigma))\n", " - mf.CovLaplace(mf.Trace(sigma))\n", " - mf.InnerProduct(mf.Ricci, sigma)\n", " )\n", "\n", " for m in range(8):\n", " t /= 2\n", " gfGt.vec.data = gfG.vec + t * gfsigma.vec\n", " mf_t = dg.RiemannianManifold(gfGt)\n", " term_t = mf_t.Scalar\n", "\n", " errold = err\n", " err = sqrt(Integrate((term_t - term_0 - t * variation) ** 2, mesh))\n", " print(\n", " f\"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}\"\n", " )\n", "\n", "\n", "print(\"\\nscalar omega_T:\")\n", "t = 1 / 8\n", "gfGt = GridFunction(gfsigma.space)\n", "errold = None\n", "err = None\n", "\n", "with TaskManager():\n", " term_0 = mf.Scalar * omega_T\n", " variation = (\n", " mf.CovDiv(mf.CovDiv(mf.S(sigma))) - mf.InnerProduct(mf.Einstein, sigma)\n", " ) * omega_T\n", "\n", " for m in range(8):\n", " t /= 2\n", " gfGt.vec.data = gfG.vec + t * gfsigma.vec\n", " mf_t = dg.RiemannianManifold(gfGt)\n", " term_t = mf_t.Scalar * mf_t.VolumeForm(VOL)\n", " errold = err\n", " err = sqrt(Integrate((term_t - term_0 - t * variation) ** 2, mesh))\n", " print(\n", " f\"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}\"\n", " )" ] }, { "cell_type": "markdown", "id": "09fede70", "metadata": {}, "source": [ "The mean curvature $H$ at a facet $F$ generalizes the geodesic curvature to higher dimensions and has the linearization\n", "$$\n", "\\begin{align*}\n", "\\frac{d}{dt}H &= 0.5\\,\\big(g(\\mathrm{div}_g(\\mathbb{S}_g\\sigma),n_g)-g_F(I\\!I,\\sigma|_F)+ \\mathrm{div}^F_g(\\sigma_{n_g})+H\\,\\sigma(n_g,n_g)\\big),\\\\\n", "\\frac{d}{dt}(H\\,\\omega_F) &= 0.5\\,\\big(g(\\mathrm{div}_g(\\mathbb{S}_g\\sigma),n_g)-g_F(\\mathbb{S}_{F,g}I\\!I,\\sigma)+ \\mathrm{div}^F_g(\\sigma_{n_g})+H\\,\\sigma(n_g,n_g)\\big)\\,\\omega_F,\n", "\\end{align*}\n", "$$\n", "where $g_F$ is the metric restricted to $F$, $\\mathbb{S}_{F,g}\\sigma=\\sigma_F-\\mathrm{tr}(\\sigma|_F)g_F$, $I\\,I$ the second fundamental form, and $\\mathrm{div}^F_g$ the covariant surface divergence." ] }, { "cell_type": "code", "execution_count": null, "id": "f99c8143", "metadata": {}, "outputs": [], "source": [ "print(\"\\nmean curvature:\")\n", "t = 1 / 8\n", "gfGt = GridFunction(gfsigma.space)\n", "errold = None\n", "err = None\n", "\n", "tv = specialcf.tangential(mesh.dim)\n", "nv = specialcf.normal(mesh.dim)\n", "\n", "nv_g = dg.VectorField(mf.normal)\n", "tv_g = mf.tangent\n", "\n", "with TaskManager():\n", " term_0 = mf.MeanCurvature\n", "\n", " variation = 0.5 * (\n", " -mf.InnerProduct(mf.SFF, sigma, BND)\n", " + mf.InnerProduct(mf.CovDiv(mf.S(sigma)), nv_g)\n", " + mf.CovDiv(mf.Contraction(sigma, nv_g, 0), BND)\n", " + mf.MeanCurvature * mf.InnerProduct(sigma, dg.TensorProduct(nv_g, nv_g))\n", " )\n", "\n", " for m in range(8):\n", " t /= 2\n", " gfGt.vec.data = gfG.vec + t * gfsigma.vec\n", " mf_t = dg.RiemannianManifold(gfGt)\n", " term_t = mf_t.MeanCurvature\n", " errold = err\n", " err = sqrt(\n", " Integrate(\n", " (term_t - term_0 - t * variation) ** 2 * dx(element_boundary=True), mesh\n", " )\n", " )\n", " print(\n", " f\"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}\"\n", " )\n", "\n", "\n", "print(\"\\nmean curvature omega_F:\")\n", "t = 1 / 8\n", "gfGt = GridFunction(gfsigma.space)\n", "errold = None\n", "err = None\n", "\n", "with TaskManager():\n", " term_0 = mf.MeanCurvature * omega_F\n", "\n", " variation = (\n", " 0.5\n", " * (\n", " -mf.InnerProduct(mf.S(mf.SFF, BND), sigma, BND)\n", " + mf.InnerProduct(mf.CovDiv(mf.S(sigma)), nv_g)\n", " + mf.CovDiv(mf.Contraction(sigma, nv_g, 0), BND)\n", " + mf.MeanCurvature * mf.InnerProduct(sigma, dg.TensorProduct(nv_g, nv_g))\n", " )\n", " * omega_F\n", " )\n", "\n", " for m in range(8):\n", " t /= 2\n", " gfGt.vec.data = gfG.vec + t * gfsigma.vec\n", " mf_t = dg.RiemannianManifold(gfGt)\n", " term_t = mf_t.MeanCurvature * mf_t.VolumeForm(BND)\n", " errold = err\n", " err = sqrt(\n", " Integrate(\n", " (term_t - term_0 - t * variation) ** 2 * dx(element_boundary=True), mesh\n", " )\n", " )\n", " print(\n", " f\"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}\"\n", " )" ] }, { "cell_type": "markdown", "id": "9d1c9dc2", "metadata": {}, "source": [ "The linearization of the Ricci curvature is given by\n", "$$\n", "\\frac{d}{dt}\\mathrm{Ric} = -0.5\\big(\\Delta_L\\sigma+\\nabla_g^2(\\mathrm{tr}\\sigma)-2\\,\\mathrm{def}_g(\\mathrm{div}_g\\sigma)\\big),\n", "$$\n", "where $\\nabla_g^2$ denotes the covariant Hessian, $\\mathrm{def}_g = 0.5(\\nabla_g+\\nabla_g^T$ the covariant defect operator (symmetric covariant derivative) for vector-fields or 1-forms, and $\\Delta_L\\sigma$ is the Licherowitz Laplacian defined by\n", "$$\n", "\\Delta_L\\sigma = \\Delta_g\\sigma - 2 \\big(\\mathfrak{R}_{ikjl}g^{la}\\sigma_{ab}g^{bk}dx^i\\otimes dx^j-\\mathrm{sym}(\\mathrm{Ric}_{ia}g^{ab}\\sigma_{bj}dx^i\\otimes dx^j)\\big).\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "id": "2f082f89", "metadata": {}, "outputs": [], "source": [ "print(\"Ricci curvature:\")\n", "t = 1 / 8\n", "gfGt = GridFunction(gfsigma.space)\n", "errold = None\n", "err = None\n", "\n", "with TaskManager():\n", " term_0 = mf.Ricci\n", " variation = -0.5 * (\n", " mf.LichnerowiczLaplacian(sigma)\n", " + mf.CovHesse(mf.Trace(sigma))\n", " - 2 * mf.CovDef(mf.CovDiv(sigma))\n", " )\n", "\n", " for m in range(8):\n", " t /= 2\n", " gfGt.vec.data = gfG.vec + t * gfsigma.vec\n", " mf_t = dg.RiemannianManifold(gfGt)\n", " term_t = mf_t.Ricci\n", "\n", " errold = err\n", " err = sqrt(\n", " Integrate(\n", " InnerProduct(\n", " term_t - term_0 - t * variation, term_t - term_0 - t * variation\n", " ),\n", " mesh,\n", " )\n", " )\n", " print(\n", " f\"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}\"\n", " )" ] }, { "cell_type": "markdown", "id": "7a641156", "metadata": {}, "source": [ "The linearization of the Einstein tensor $G=J\\mathrm{Ric}:=\\mathrm{Ric}-0.5\\,S\\,g$ follows by the product rule and simplifying\n", "$$\n", "\\begin{align*}\n", "\\frac{d}{dt}G &= -0.5\\Big(\\Delta_L\\sigma - 2\\mathrm{def}_g(\\mathrm{div}_g\\sigma)+\\big(\\mathrm{div}_g\\mathrm{div}_g(J(\\sigma))-g(\\mathrm{Ric},\\sigma)\\big)\\,g + S\\,\\sigma\\Big),\\\\\n", "\\frac{d}{dt}(G\\,\\omega_T) &= 0.5\\Big(2\\mathrm{ein}_g\\sigma +2\\mathfrak{R}_{ikjl}g^{la}\\sigma_{ab}g^{bk} + g(\\mathrm{Ric},\\sigma)\\,g+\\mathrm{tr}\\sigma\\,\\mathrm{Ric}+\\mathrm{sym}(\\mathrm{Ric}_{ia}g^{ab}\\sigma_{bj})-S\\,\\sigma-0.5\\mathrm{tr}(\\sigma)\\,S\\,g \\Big)\\,\\omega_T.\n", "\\end{align*}\n", "$$\n", "Here, $\\mathrm{ein}_g\\sigma= J\\mathrm{def}_g(\\mathrm{div}_gJ\\sigma)-0.5\\Delta_gJ\\sigma$ is the covariant linearized Einstein operator." ] }, { "cell_type": "markdown", "id": "07add940", "metadata": {}, "source": [ "##" ] }, { "cell_type": "code", "execution_count": null, "id": "e87638c1", "metadata": {}, "outputs": [], "source": [ "print(\"Einstein tensor:\")\n", "t = 1 / 8\n", "gfGt = GridFunction(gfsigma.space)\n", "errold = None\n", "err = None\n", "\n", "with TaskManager():\n", " term_0 = mf.Einstein\n", " variation = -0.5 * (\n", " mf.LichnerowiczLaplacian(sigma)\n", " - 2 * mf.CovDef(mf.CovDiv(mf.J(sigma)))\n", " + (mf.CovDiv(mf.CovDiv(mf.S(sigma))) - mf.InnerProduct(mf.Ricci, sigma)) * gfG\n", " + mf.Scalar * sigma\n", " )\n", "\n", " for m in range(6):\n", " t /= 2\n", " gfGt.vec.data = gfG.vec + t * gfsigma.vec\n", " mf_t = dg.RiemannianManifold(gfGt)\n", " term_t = mf_t.Einstein\n", "\n", " errold = err\n", " err = sqrt(\n", " Integrate(\n", " InnerProduct(\n", " term_t - term_0 - t * variation, term_t - term_0 - t * variation\n", " ),\n", " mesh,\n", " )\n", " )\n", " print(\n", " f\"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}\"\n", " )\n", "\n", "\n", "print(\"\\nEinstein omega_T:\")\n", "t = 1 / 8\n", "gfGt = GridFunction(gfsigma.space)\n", "errold = None\n", "err = None\n", "\n", "with TaskManager():\n", " term_0 = mf.Einstein * omega_T\n", " variation = (\n", " mf.CovEin(sigma)\n", " + fem.Einsum(\"ikjl,lk->ij\", mf.Riemann, Inv(gfG) * sigma * Inv(gfG))\n", " + 0.5 * mf.InnerProduct(mf.Ricci, sigma) * gfG\n", " + 0.5 * mf.Trace(sigma) * mf.Ricci\n", " + Sym(mf.Ricci * Inv(gfG) * sigma)\n", " - 0.5 * mf.Scalar * sigma\n", " - 0.25 * mf.Trace(sigma) * mf.Scalar * gfG\n", " ) * omega_T\n", "\n", " for m in range(6):\n", " t /= 2\n", " gfGt.vec.data = gfG.vec + t * gfsigma.vec\n", " mf_t = dg.RiemannianManifold(gfGt)\n", " term_t = mf_t.Einstein * mf_t.VolumeForm(VOL)\n", "\n", " errold = err\n", " err = sqrt(\n", " Integrate(\n", " InnerProduct(\n", " term_t - term_0 - t * variation, term_t - term_0 - t * variation\n", " ),\n", " mesh,\n", " )\n", " )\n", " print(\n", " f\"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}\"\n", " )\n" ] }, { "cell_type": "markdown", "id": "476d3bf5", "metadata": {}, "source": [ "The linearization of the full Riemann curvature tensor is given by\n", "$$\n", "\\frac{d}{dt}\\mathfrak{R} = -2\\,\\mathrm{Inc}_g\\sigma +0.5\\Big(\\mathfrak{R}_{ijal}g^{ab}\\sigma_{bk}+\\mathfrak{R}_{ijka}g^{ab}\\sigma_{bl}\\Big)dx^i\\otimes dx^j\\otimes dx^k\\otimes dx^l,\n", "$$\n", "where \n", "$$\n", "(\\mathrm{Inc}_g\\sigma)(X,Y,Z,W) = -\\frac{1}{4}\\big(\\nabla^2_{X,Z}\\sigma(Y,W)-\\nabla^2_{Y,Z}\\sigma(X,W)-\\nabla^2_{X,W}\\sigma(Y,Z)+\\nabla^2_{Y,W}\\sigma(X,Z)\\big)\n", "$$\n", "is covariant incompatibility operator in arbitrary dimensions." ] }, { "cell_type": "code", "execution_count": null, "id": "21c8ccf2", "metadata": {}, "outputs": [], "source": [ "print(\"\\nRiemann curvature tensor:\")\n", "t = 1 / 8\n", "gfGt = GridFunction(gfsigma.space)\n", "errold = None\n", "err = None\n", "\n", "with TaskManager():\n", " term_0 = mf.Riemann\n", " variation = -2 * mf.CovInc(sigma) + 0.5 * (\n", " fem.Einsum(\"ijal,ak->ijkl\", mf.Riemann, Inv(gfG) * sigma)\n", " + fem.Einsum(\"ijka,al->ijkl\", mf.Riemann, Inv(gfG) * sigma)\n", " )\n", "\n", " for m in range(8):\n", " t /= 2\n", " gfGt.vec.data = gfG.vec + t * gfsigma.vec\n", " mf_t = dg.RiemannianManifold(gfGt)\n", " term_t = mf_t.Riemann\n", "\n", " errold = err\n", " err = sqrt(\n", " Integrate(\n", " InnerProduct(\n", " term_t - term_0 - t * variation, term_t - term_0 - t * variation\n", " ),\n", " mesh,\n", " )\n", " )\n", " print(\n", " f\"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}\"\n", " )" ] }, { "cell_type": "markdown", "id": "1a008f27", "metadata": {}, "source": [ "The second fundamental form $I\\!I(X,Y) = -g(\\nabla_X\\nu,Y)$, $X,Y\\in \\mathfrak{X}(F)$ at facet $F$ has the following linearization\n", "$$\n", "\\Big(\\frac{d}{dt}I\\!I\\Big)(X,Y) = 0.5\\,\\sigma(n_g,n_g)I\\!I(X,Y)+0.5\\big((\\nabla_X\\sigma)(n_g,Y)+(\\nabla_Y\\sigma)(n_g,X)-(\\nabla_{n_g}\\sigma)(X,Y)\\big),\\qquad X,Y\\in\\mathfrak{X}(F)\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "id": "2430f9c7", "metadata": {}, "outputs": [], "source": [ "print(\"second fundamental form:\")\n", "t = 1 / 8\n", "gfGt = GridFunction(gfsigma.space)\n", "errold = None\n", "err = None\n", "\n", "Q_eucl = Id(mesh.dim) - OuterProduct(nv, nv)\n", "\n", "with TaskManager():\n", " term_0 = Q_eucl * mf.SFF * Q_eucl\n", " variation = (\n", " Q_eucl\n", " * (\n", " 0.5 * mf.InnerProduct(sigma, dg.TensorProduct(nv_g, nv_g)) * mf.SFF\n", " + 0.5\n", " * (\n", " 2 * Sym(mf.Contraction(mf.CovDeriv(sigma), nv_g, 1))\n", " - mf.Contraction(mf.CovDeriv(sigma), nv_g, 0)\n", " )\n", " )\n", " * Q_eucl\n", " )\n", "\n", " for m in range(8):\n", " t /= 2\n", " gfGt.vec.data = gfG.vec + t * gfsigma.vec\n", " mf_t = dg.RiemannianManifold(gfGt)\n", " term_t = Q_eucl * mf_t.SFF * Q_eucl\n", "\n", " errold = err\n", " err = sqrt(\n", " Integrate(\n", " InnerProduct(\n", " term_t - term_0 - t * variation, term_t - term_0 - t * variation\n", " )\n", " * dx(element_boundary=True),\n", " mesh,\n", " )\n", " )\n", " print(\n", " f\"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}\"\n", " )" ] }, { "cell_type": "code", "execution_count": null, "id": "13deb195", "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 }