{ "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, normal_sign=-1)\n", "\n", "# quadratic convergence rate is expected\n", "rate_tol = 1.75" ] }, { "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, normal_sign=-1)\n", " term_t = mf_t.VolumeForm(VOL)\n", "\n", " errold = err\n", " err = sqrt(Integrate((term_t - term_0 - t * variation) ** 2, mesh))\n", " rate = round(log2(errold / err), 3) if errold else \"-\"\n", " print(f\"err = {err:.10f}, order = {rate}\")\n", " assert rate == \"-\" or rate > rate_tol\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, normal_sign=-1)\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", " rate = round(log2(errold / err), 3) if errold else \"-\"\n", " print(f\"err = {err:.10f}, order = {rate}\")\n", " assert rate == \"-\" or rate > rate_tol" ] }, { "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, normal_sign=-1)\n", " term_t = mf_t.Gauss\n", "\n", " errold = err\n", " err = sqrt(Integrate((term_t - term_0 - t * variation) ** 2, mesh))\n", " rate = round(log2(errold / err), 3) if errold else \"-\"\n", " print(f\"err = {err:.10f}, order = {rate}\")\n", " assert rate == \"-\" or rate > rate_tol\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, normal_sign=-1)\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", " rate = round(log2(errold / err), 3) if errold else \"-\"\n", " print(f\"err = {err:.10f}, order = {rate}\")\n", " assert rate == \"-\" or rate > rate_tol" ] }, { "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 = 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(mf.InnerProduct(sigma, dg.TensorProduct(nv_g, tv_g))),\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, normal_sign=-1)\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", " rate = round(log2(errold / err), 3) if errold else \"-\"\n", " print(f\"err = {err:.10f}, order = {rate}\")\n", " assert rate == \"-\" or rate > rate_tol" ] }, { "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, normal_sign=-1)" ] }, { "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, normal_sign=-1)\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", " rate = round(log2(errold / err), 3) if errold else \"-\"\n", " print(f\"err = {err:.10f}, order = {rate}\")\n", " assert rate == \"-\" or rate > rate_tol" ] }, { "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", "$$\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", "\n", "where $\\Delta_g$ is the Laplace-Beltrami operator, $\\mathrm{Ric}$ the Ricci curvature tensor, and $G$ the Einstein tensor.\n" ] }, { "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, normal_sign=-1)\n", " term_t = mf_t.Scalar\n", "\n", " errold = err\n", " err = sqrt(Integrate((term_t - term_0 - t * variation) ** 2, mesh))\n", " rate = round(log2(errold / err), 3) if errold else \"-\"\n", " print(f\"err = {err:.10f}, order = {rate}\")\n", " assert rate == \"-\" or rate > rate_tol\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", " rate = round(log2(errold / err), 3) if errold else \"-\"\n", " print(f\"err = {err:.10f}, order = {rate}\")\n", " assert rate == \"-\" or rate > rate_tol" ] }, { "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", "$$\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", "\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.\n" ] }, { "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 = 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, normal_sign=-1)\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", " rate = round(log2(errold / err), 3) if errold else \"-\"\n", " print(f\"err = {err:.10f}, order = {rate}\")\n", " assert rate == \"-\" or rate > rate_tol\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, normal_sign=-1)\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", " rate = round(log2(errold / err), 3) if errold else \"-\"\n", " print(f\"err = {err:.10f}, order = {rate}\")\n", " assert rate == \"-\" or rate > rate_tol" ] }, { "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)- 2\\mathfrak{R}_{aijb}\\sigma^{ab}dx^i\\otimes dx^j,\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}\\sigma^{lk}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", "\n", " variation = -0.5 * (\n", " mf.LichnerowiczLaplacian(sigma)\n", " + mf.CovHesse(mf.Trace(sigma))\n", " - 2 * mf.CovDef(mf.CovDiv(sigma))\n", " ) - 2 * dg.TensorField(\n", " fem.Einsum(\"ijkl,il->jk\", mf.Riemann, Inv(gfG) * sigma * Inv(gfG)),\n", " \"11\",\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, normal_sign=-1)\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", " rate = round(log2(errold / err), 3) if errold else \"-\"\n", " print(f\"err = {err:.10f}, order = {rate}\")\n", " assert rate == \"-\" or rate > rate_tol" ] }, { "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", "$$\n", "\\begin{align*}\n", "\\frac{d}{dt}G &= -0.5\\Big(\\Delta_L\\sigma + 4\\mathfrak{R}_{aijb}\\sigma^{ab}dx^i\\otimes dx^j - 2\\mathrm{def}_g(\\mathrm{div}_gJ(\\sigma))+\\big(\\mathrm{div}_g\\mathrm{div}_g(S(\\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", "\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.\n" ] }, { "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", "\n", " variation = -0.5 * (\n", " mf.LichnerowiczLaplacian(sigma)\n", " + 4\n", " * dg.TensorField(\n", " fem.Einsum(\"ijkl,il->jk\", mf.Riemann, Inv(gfG) * sigma * Inv(gfG)),\n", " \"11\",\n", " )\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(8):\n", " t /= 2\n", " gfGt.vec.data = gfG.vec + t * gfsigma.vec\n", " mf_t = dg.RiemannianManifold(gfGt, normal_sign=-1)\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", " rate = round(log2(errold / err), 3) if errold else \"-\"\n", " print(f\"err = {err:.10f}, order = {rate}\")\n", " assert rate == \"-\" or rate > rate_tol\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.coef * 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(8):\n", " t /= 2\n", " gfGt.vec.data = gfG.vec + t * gfsigma.vec\n", " mf_t = dg.RiemannianManifold(gfGt, normal_sign=-1)\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", " rate = round(log2(errold / err), 3) if errold else \"-\"\n", " print(f\"err = {err:.10f}, order = {rate}\")\n", " assert rate == \"-\" or rate > rate_tol" ] }, { "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, normal_sign=-1)\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", " rate = round(log2(errold / err), 3) if errold else \"-\"\n", " print(f\"err = {err:.10f}, order = {rate}\")\n", " assert rate == \"-\" or rate > rate_tol" ] }, { "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", "$$\n", "Note, that we project it to the tangent space of $F$ using the Euclidean normal vector for comparison, as $I\\!I$ only acts on the tangent space." ] }, { "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(specialcf.normal(3), specialcf.normal(3))\n", "nv_g = mf.normal\n", "\n", "with TaskManager():\n", " term_0 = Q_eucl * mf.SFF.coef * Q_eucl\n", "\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(6):\n", " t /= 2\n", " gfGt.vec.data = gfG.vec + t * gfsigma.vec\n", " mf_t = dg.RiemannianManifold(gfGt, normal_sign=-1)\n", " term_t = Q_eucl * mf_t.SFF.coef * 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", " rate = round(log2(errold / err), 3) if errold else \"-\"\n", " print(f\"err = {err:.10f}, order = {rate}\")\n", " assert rate == \"-\" or rate > rate_tol" ] }, { "cell_type": "code", "execution_count": null, "id": "13deb195", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "ngsolve-dev (3.12.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.12.3" } }, "nbformat": 4, "nbformat_minor": 5 }