{ "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": [ "### Volume forms\n", "\n", "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": [ "### Gauss curvature\n", "\n", "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": [ "### Geodesic curvature\n", "\n", "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.CovDerivative(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", "sigma_form = dg.DoubleForm(sigma, p=1, q=1, dim=mesh.dim)\n", "\n", "mf = dg.RiemannianManifold(gfG, normal_sign=-1)" ] }, { "cell_type": "markdown", "id": "489e2d11", "metadata": {}, "source": [ "### Edge volume form\n", "\n", "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": [ "### Scalar curvature\n", "\n", "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": [ "### Mean curvature\n", "\n", "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": [ "### Ricci curvature\n", "\n", "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 deformation operator (symmetric covariant derivative) for vector fields or 1-forms, and $\\Delta_L\\sigma$ is the Lichnerowicz 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.CovHessian(mf.Trace(sigma))\n", " - 2 * mf.CovDef(mf.CovDiv(sigma))\n", " ) + 2 * mf.TraceSigma(mf.Riemann, sigma_form)\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": [ "### Einstein tensor\n", "\n", "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": "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 * mf.TraceSigma(mf.Riemann, sigma_form)\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", " + mf.TraceSigma(mf.Riemann, sigma_form)\n", " + 0.5 * mf.InnerProduct(mf.Ricci, sigma) * gfG\n", " + 0.5 * mf.Trace(sigma) * mf.Ricci\n", " + Sym(mf.Ricci.coef * mf.Raise(sigma, 0).coef)\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": "07add940", "metadata": {}, "source": [ "### Riemann curvature tensor" ] }, { "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, mf.Raise(sigma, 0).coef)\n", " + fem.Einsum(\"ijka,al->ijkl\", mf.Riemann, mf.Raise(sigma, 0).coef)\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": [ "### Second fundamental form\n", "\n", "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", "The formula holds for tangent vectors to the fixed facet $F$. For comparison, we restrict all three tensors using the same Euclidean reference manifold: `euclidean.ProjectTensor(..., mode=\"F\")` applies the fixed Euclidean tangential projector to both arguments. This projector is independent of the varied metric $g+t\\sigma$." ] }, { "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", "euclidean = dg.RiemannianManifold(Id(mesh.dim))\n", "nv_g = mf.normal\n", "\n", "with TaskManager():\n", " term_0 = euclidean.ProjectTensor(mf.SFF, mode=\"F\")\n", "\n", " cov_sigma = mf.CovDerivative(sigma)\n", " A = mf.Contraction(cov_sigma, nv_g, 1)\n", " B = mf.Contraction(cov_sigma, nv_g, 0)\n", " sigma_nn = mf.InnerProduct(sigma, dg.TensorProduct(nv_g, nv_g))\n", " variation = euclidean.ProjectTensor(\n", " 0.5 * sigma_nn * mf.SFF + 0.5 * (2 * dg.Sym(A) - B),\n", " mode=\"F\",\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 = euclidean.ProjectTensor(mf_t.SFF, mode=\"F\")\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 }