Linearization of curvature quantities
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}\)
\[F(g+t\sigma)-F(g)-t\, \dot{F}(g),\]
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\).
2D
[1]:
from math import log2
from ngsolve import *
from netgen.occ import *
import ngsdiffgeo as dg
from ngsolve.meshes import MakeStructured2DMesh
mesh = MakeStructured2DMesh(False, nx=6, ny=6)
peak = 0.5 * x**2 - 1 / 12 * x**4 + 0.5 * y**2 - 1 / 12 * y**4
F = CF((1, 0, 0, 1, peak.Diff(x), peak.Diff(y)), dims=(3, 2))
Gex = F.trans * F
cfsigma = dg.TensorField(
Sym(
CF(
(
1 + 10 * x * y**3 - x**2,
0.2 * y**4 * x - y,
0.2 * sin(x * y),
1 + cos(x) * y**2,
),
dims=(2, 2),
)
),
"11",
)
order = 3
gfG = GridFunction(HCurlCurl(mesh, order=order))
gfsigma = GridFunction(HCurlCurl(mesh, order=order))
with TaskManager():
gfG.Set(Gex)
gfsigma.Set(cfsigma)
sigma = dg.TensorField(gfsigma, "11")
mf = dg.RiemannianManifold(gfG)
We start by verifying the variations of the volume forms
\[\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\]
[2]:
t = 1 / 8
gfGt = GridFunction(gfsigma.space)
errold = None
err = None
omega_T = mf.VolumeForm(VOL)
omega_F = mf.VolumeForm(BND)
print("omega_T:")
with TaskManager():
term_0 = omega_T
variation = 0.5 * mf.Trace(sigma) * omega_T
for m in range(8):
t /= 2
gfGt.vec.data = gfG.vec + t * gfsigma.vec
mf_t = dg.RiemannianManifold(gfGt)
term_t = mf_t.VolumeForm(VOL)
errold = err
err = sqrt(Integrate((term_t - term_0 - t * variation) ** 2, mesh))
print(
f"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}"
)
print("\nomega_F:")
t = 1 / 8
errold = None
err = None
with TaskManager():
term_0 = omega_F
variation = 0.5 * mf.Trace(sigma, BND) * omega_F
for m in range(8):
t /= 2
gfGt.vec.data = gfG.vec + t * gfsigma.vec
mf_t = dg.RiemannianManifold(gfGt)
term_t = mf_t.VolumeForm(BND)
errold = err
err = sqrt(
Integrate(
(term_t - term_0 - t * variation) ** 2 * dx(element_boundary=True), mesh
)
)
print(
f"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}"
)
omega_T:
err = 0.0035249629, order = -
err = 0.0009621393, order = 1.873
err = 0.0002523502, order = 1.931
err = 0.0000646952, order = 1.964
err = 0.0000163839, order = 1.981
err = 0.0000041228, order = 1.991
err = 0.0000010341, order = 1.995
err = 0.0000002590, order = 1.998
omega_F:
err = 0.0189513367, order = -
err = 0.0050394053, order = 1.911
err = 0.0013028271, order = 1.952
err = 0.0003314756, order = 1.975
err = 0.0000836174, order = 1.987
err = 0.0000209997, order = 1.993
err = 0.0000052620, order = 1.997
err = 0.0000013170, order = 1.998
Next, we look at the Gauss curvature \(K\) and the Gauss curvature multiplied by the volume form \(K\,\omega_T\). The variations are
\[\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,\]
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\).
[3]:
print("Gauss curvature:")
t = 1 / 8
gfGt = GridFunction(gfsigma.space)
errold = None
err = None
with TaskManager():
term_0 = mf.Gauss
variation = 0.5 * (
mf.CovDiv(mf.CovDiv(sigma))
- mf.CovLaplace(mf.Trace(sigma))
- mf.InnerProduct(sigma, mf.Ricci)
)
for m in range(8):
t /= 2
gfGt.vec.data = gfG.vec + t * gfsigma.vec
mf_t = dg.RiemannianManifold(gfGt)
term_t = mf_t.Gauss
errold = err
err = sqrt(Integrate((term_t - term_0 - t * variation) ** 2, mesh))
print(
f"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}"
)
print("\nK omega_T:")
t = 1 / 8
gfGt = GridFunction(gfsigma.space)
errold = None
err = None
with TaskManager():
term_0 = mf.Gauss * omega_T
variation = 0.5 * mf.CovDiv(mf.CovDiv(mf.S(sigma))) * omega_T
for m in range(8):
t /= 2
gfGt.vec.data = gfG.vec + t * gfsigma.vec
mf_t = dg.RiemannianManifold(gfGt)
term_t = mf_t.Gauss * mf_t.VolumeForm(VOL)
errold = err
err = sqrt(Integrate((term_t - term_0 - t * variation) ** 2, mesh))
print(
f"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}"
)
Gauss curvature:
err = 0.1145910933, order = -
err = 0.0339083919, order = 1.757
err = 0.0093398657, order = 1.86
err = 0.0024606599, order = 1.924
err = 0.0006322177, order = 1.961
err = 0.0001602784, order = 1.98
err = 0.0000403537, order = 1.99
err = 0.0000101243, order = 1.995
K omega_T:
err = 0.1108513348, order = -
err = 0.0320222891, order = 1.791
err = 0.0086919313, order = 1.881
err = 0.0022712561, order = 1.936
err = 0.0005810201, order = 1.967
err = 0.0001469688, order = 1.983
err = 0.0000369605, order = 1.991
err = 0.0000092677, order = 1.996
For the geodesic curvature \(\kappa\) at edges we have that
\[\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,\]
where \(n_g\) and \(t_g\) are the \(g\)-normal and tangential vectors.
[4]:
print("kappa omega_F:")
t = 1 / 8
gfGt = GridFunction(gfsigma.space)
errold = None
err = None
tv = specialcf.tangential(mesh.dim)
nv = specialcf.normal(mesh.dim)
nv_g = dg.VectorField(mf.normal)
tv_g = mf.tangent
with TaskManager():
term_0 = mf.GeodesicCurvature * omega_F
variation = (
0.5
* (
mf.InnerProduct(mf.CovDiv(mf.S(sigma)), nv_g)
+ mf.InnerProduct(
mf.CovDeriv(
dg.ScalarField(mf.InnerProduct(sigma, dg.TensorProduct(nv_g, tv_g)))
),
tv_g,
)
)
* omega_F
)
for m in range(8):
t /= 2
gfGt.vec.data = gfG.vec + t * gfsigma.vec
mf_t = dg.RiemannianManifold(gfGt)
term_t = mf_t.GeodesicCurvature * mf_t.VolumeForm(BND)
errold = err
err = sqrt(
Integrate(
(term_t - term_0 - t * variation) ** 2 * dx(element_boundary=True), mesh
)
)
print(
f"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}"
)
kappa omega_F:
err = 0.1681449576, order = -
err = 0.0469285654, order = 1.841
err = 0.0124799881, order = 1.911
err = 0.0032245567, order = 1.952
err = 0.0008200097, order = 1.975
err = 0.0002067903, order = 1.987
err = 0.0000519246, order = 1.994
err = 0.0000130097, order = 1.997
3D
[5]:
from math import log2
from ngsolve import *
from ngsolve.fem import Einsum
from netgen.occ import *
import ngsdiffgeo as dg
from ngsolve.meshes import MakeStructured3DMesh
mesh = MakeStructured3DMesh(False, nx=2, ny=2, nz=2)
peak = (
0.5 * x**2 - 1 / 12 * x**4 + 0.5 * y**2 - 1 / 12 * y**4 + 0.5 * z**2 - 1 / 12 * z**4
)
F = CF(
(1, 0, 0, 0, 1, 0, 0, 0, 1, peak.Diff(x), peak.Diff(y), peak.Diff(z)),
dims=(4, 3),
)
Gex = F.trans * F
cfsigma = dg.TensorField(
Sym(
CF(
(
1 + 10 * x * y**3 - x**2,
0.2 * y**4 * x - y,
0.2 * z**4 * x - z,
0.2 * sin(x * y),
1 + cos(x) * y**2,
0.2 * sin(x * z),
0.2 * z**4 * x - z,
0.2 * sin(x * z),
1 + cos(x) * z**2,
),
dims=(3, 3),
)
),
"11",
)
order = 2
gfG = GridFunction(HCurlCurl(mesh, order=order))
gfsigma = GridFunction(HCurlCurl(mesh, order=order))
with TaskManager():
gfG.Set(Gex)
gfsigma.Set(cfsigma)
sigma = dg.TensorField(gfsigma, "11")
mf = dg.RiemannianManifold(gfG)
In 3D, we also have the volume form on codimension 2 boundaries
\[\frac{d}{dt}(\omega_E) = 0.5\,\mathrm{tr}(\sigma|_E)\,\omega_E.\]
[6]:
t = 1 / 8
gfGt = GridFunction(gfsigma.space)
errold = None
err = None
omega_T = mf.VolumeForm(VOL)
omega_F = mf.VolumeForm(BND)
omega_E = mf.VolumeForm(BBND)
print("omega_E:")
t = 1 / 8
errold = None
err = None
with TaskManager():
term_0 = omega_E
variation = 0.5 * mf.Trace(sigma, BBND) * omega_E
for m in range(10):
t /= 2
gfGt.vec.data = gfG.vec + t * gfsigma.vec
mf_t = dg.RiemannianManifold(gfGt)
term_t = mf_t.VolumeForm(BBND)
f = BilinearForm(gfsigma.space)
f += Variation(
InnerProduct(
term_t - term_0 - t * variation, term_t - term_0 - t * variation
)
* dx(element_vb=BBND)
)
errold = err
err = sqrt(f.Energy(gfsigma.vec))
print(
f"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}"
)
omega_E:
err = 0.0321053102, order = -
err = 0.0086069850, order = 1.899
err = 0.0022355777, order = 1.945
err = 0.0005702281, order = 1.971
err = 0.0001440331, order = 1.985
err = 0.0000361967, order = 1.992
err = 0.0000090730, order = 1.996
err = 0.0000022712, order = 1.998
err = 0.0000005682, order = 1.999
err = 0.0000001421, order = 2.0
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
\[\begin{split}\begin{align*}
\frac{d}{dt}S &= \mathrm{div}_g\mathrm{div}_g\sigma-\Delta_g\mathrm{tr}\sigma-g(\mathrm{Ric},\sigma),\\
\frac{d}{dt}(S\,\omega_T) &= \big(\mathrm{div}_g\mathrm{div}_g(\mathbb{S}_g\sigma)-g(\sigma,G)\big)\,\omega_T,
\end{align*}\end{split}\]
where \(\Delta_g\) is the Laplace-Beltrami operator, \(\mathrm{Ric}\) the Ricci curvature tensor, and \(G\) the Einstein tensor.
[7]:
print("\nscalar:")
t = 1 / 8
gfGt = GridFunction(gfsigma.space)
errold = None
err = None
with TaskManager():
term_0 = mf.Scalar
variation = (
mf.CovDiv(mf.CovDiv(sigma))
- mf.CovLaplace(mf.Trace(sigma))
- mf.InnerProduct(mf.Ricci, sigma)
)
for m in range(8):
t /= 2
gfGt.vec.data = gfG.vec + t * gfsigma.vec
mf_t = dg.RiemannianManifold(gfGt)
term_t = mf_t.Scalar
errold = err
err = sqrt(Integrate((term_t - term_0 - t * variation) ** 2, mesh))
print(
f"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}"
)
print("\nscalar omega_T:")
t = 1 / 8
gfGt = GridFunction(gfsigma.space)
errold = None
err = None
with TaskManager():
term_0 = mf.Scalar * omega_T
variation = (
mf.CovDiv(mf.CovDiv(mf.S(sigma))) - mf.InnerProduct(mf.Einstein, sigma)
) * omega_T
for m in range(8):
t /= 2
gfGt.vec.data = gfG.vec + t * gfsigma.vec
mf_t = dg.RiemannianManifold(gfGt)
term_t = mf_t.Scalar * mf_t.VolumeForm(VOL)
errold = err
err = sqrt(Integrate((term_t - term_0 - t * variation) ** 2, mesh))
print(
f"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}"
)
scalar:
err = 0.2264855748, order = -
err = 0.0662693640, order = 1.773
err = 0.0181455708, order = 1.869
err = 0.0047662191, order = 1.929
err = 0.0012227375, order = 1.963
err = 0.0003097517, order = 1.981
err = 0.0000779575, order = 1.99
err = 0.0000195550, order = 1.995
scalar omega_T:
err = 0.1965983500, order = -
err = 0.0568023431, order = 1.791
err = 0.0154271989, order = 1.88
err = 0.0040331370, order = 1.936
err = 0.0010320349, order = 1.966
err = 0.0002610946, order = 1.983
err = 0.0000656671, order = 1.991
err = 0.0000164664, order = 1.996
The mean curvature \(H\) at a facet \(F\) generalizes the geodesic curvature to higher dimensions and has the linearization
\[\begin{split}\begin{align*}
\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),\\
\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,
\end{align*}\end{split}\]
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.
[8]:
print("\nmean curvature:")
t = 1 / 8
gfGt = GridFunction(gfsigma.space)
errold = None
err = None
tv = specialcf.tangential(mesh.dim)
nv = specialcf.normal(mesh.dim)
nv_g = dg.VectorField(mf.normal)
tv_g = mf.tangent
with TaskManager():
term_0 = mf.MeanCurvature
variation = 0.5 * (
-mf.InnerProduct(mf.SFF, sigma, BND)
+ mf.InnerProduct(mf.CovDiv(mf.S(sigma)), nv_g)
+ mf.CovDiv(mf.Contraction(sigma, nv_g, 0), BND)
+ mf.MeanCurvature * mf.InnerProduct(sigma, dg.TensorProduct(nv_g, nv_g))
)
for m in range(8):
t /= 2
gfGt.vec.data = gfG.vec + t * gfsigma.vec
mf_t = dg.RiemannianManifold(gfGt)
term_t = mf_t.MeanCurvature
errold = err
err = sqrt(
Integrate(
(term_t - term_0 - t * variation) ** 2 * dx(element_boundary=True), mesh
)
)
print(
f"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}"
)
print("\nmean curvature omega_F:")
t = 1 / 8
gfGt = GridFunction(gfsigma.space)
errold = None
err = None
with TaskManager():
term_0 = mf.MeanCurvature * omega_F
variation = (
0.5
* (
-mf.InnerProduct(mf.S(mf.SFF, BND), sigma, BND)
+ mf.InnerProduct(mf.CovDiv(mf.S(sigma)), nv_g)
+ mf.CovDiv(mf.Contraction(sigma, nv_g, 0), BND)
+ mf.MeanCurvature * mf.InnerProduct(sigma, dg.TensorProduct(nv_g, nv_g))
)
* omega_F
)
for m in range(8):
t /= 2
gfGt.vec.data = gfG.vec + t * gfsigma.vec
mf_t = dg.RiemannianManifold(gfGt)
term_t = mf_t.MeanCurvature * mf_t.VolumeForm(BND)
errold = err
err = sqrt(
Integrate(
(term_t - term_0 - t * variation) ** 2 * dx(element_boundary=True), mesh
)
)
print(
f"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}"
)
mean curvature:
err = 0.1932558440, order = -
err = 0.0745551239, order = 1.374
err = 0.0335152669, order = 1.153
err = 0.0164625840, order = 1.026
err = 0.0082781158, order = 0.992
err = 0.0041682973, order = 0.99
err = 0.0020937872, order = 0.993
err = 0.0010496021, order = 0.996
mean curvature omega_F:
err = 0.1905932278, order = -
err = 0.0819944202, order = 1.217
err = 0.0390041283, order = 1.072
err = 0.0193457397, order = 1.012
err = 0.0096915319, order = 0.997
err = 0.0048586616, order = 0.996
err = 0.0024336473, order = 0.997
err = 0.0012180433, order = 0.999
The linearization of the Ricci curvature is given by
\[\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),\]
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
\[\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).\]
[9]:
print("Ricci curvature:")
t = 1 / 8
gfGt = GridFunction(gfsigma.space)
errold = None
err = None
with TaskManager():
term_0 = mf.Ricci
variation = -0.5 * (
mf.LichnerowiczLaplacian(sigma)
+ mf.CovHesse(mf.Trace(sigma))
- 2 * mf.CovDef(mf.CovDiv(sigma))
)
for m in range(8):
t /= 2
gfGt.vec.data = gfG.vec + t * gfsigma.vec
mf_t = dg.RiemannianManifold(gfGt)
term_t = mf_t.Ricci
errold = err
err = sqrt(
Integrate(
InnerProduct(
term_t - term_0 - t * variation, term_t - term_0 - t * variation
),
mesh,
)
)
print(
f"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}"
)
Ricci curvature:
err = 0.1749835678, order = -
err = 0.0511596234, order = 1.774
err = 0.0139985277, order = 1.87
err = 0.0036752126, order = 1.929
err = 0.0009425917, order = 1.963
err = 0.0002387484, order = 1.981
err = 0.0000600830, order = 1.99
err = 0.0000150708, order = 1.995
The linearization of the Einstein tensor \(G=J\mathrm{Ric}:=\mathrm{Ric}-0.5\,S\,g\) follows by the product rule and simplifying
\[\begin{split}\begin{align*}
\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),\\
\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.
\end{align*}\end{split}\]
Here, \(\mathrm{ein}_g\sigma= J\mathrm{def}_g(\mathrm{div}_gJ\sigma)-0.5\Delta_gJ\sigma\) is the covariant linearized Einstein operator.
[10]:
print("Einstein tensor:")
t = 1 / 8
gfGt = GridFunction(gfsigma.space)
errold = None
err = None
with TaskManager():
term_0 = mf.Einstein
variation = -0.5 * (
mf.LichnerowiczLaplacian(sigma)
- 2 * mf.CovDef(mf.CovDiv(mf.J(sigma)))
+ (mf.CovDiv(mf.CovDiv(mf.S(sigma))) - mf.InnerProduct(mf.Ricci, sigma)) * gfG
+ mf.Scalar * sigma
)
for m in range(6):
t /= 2
gfGt.vec.data = gfG.vec + t * gfsigma.vec
mf_t = dg.RiemannianManifold(gfGt)
term_t = mf_t.Einstein
errold = err
err = sqrt(
Integrate(
InnerProduct(
term_t - term_0 - t * variation, term_t - term_0 - t * variation
),
mesh,
)
)
print(
f"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}"
)
print("\nEinstein omega_T:")
t = 1 / 8
gfGt = GridFunction(gfsigma.space)
errold = None
err = None
with TaskManager():
term_0 = mf.Einstein * omega_T
variation = (
mf.CovEin(sigma)
+ fem.Einsum("ikjl,lk->ij", mf.Riemann, Inv(gfG) * sigma * Inv(gfG))
+ 0.5 * mf.InnerProduct(mf.Ricci, sigma) * gfG
+ 0.5 * mf.Trace(sigma) * mf.Ricci
+ Sym(mf.Ricci * Inv(gfG) * sigma)
- 0.5 * mf.Scalar * sigma
- 0.25 * mf.Trace(sigma) * mf.Scalar * gfG
) * omega_T
for m in range(6):
t /= 2
gfGt.vec.data = gfG.vec + t * gfsigma.vec
mf_t = dg.RiemannianManifold(gfGt)
term_t = mf_t.Einstein * mf_t.VolumeForm(VOL)
errold = err
err = sqrt(
Integrate(
InnerProduct(
term_t - term_0 - t * variation, term_t - term_0 - t * variation
),
mesh,
)
)
print(
f"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}"
)
Einstein tensor:
err = 0.1512166526, order = -
err = 0.0443818474, order = 1.769
err = 0.0121730979, order = 1.866
err = 0.0032003156, order = 1.927
err = 0.0008213934, order = 1.962
err = 0.0002081289, order = 1.981
Einstein omega_T:
err = 0.1340645474, order = -
err = 0.0386422341, order = 1.795
err = 0.0104839120, order = 1.882
err = 0.0027394745, order = 1.936
err = 0.0007008378, order = 1.967
err = 0.0001772848, order = 1.983
The linearization of the full Riemann curvature tensor is given by
\[\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,\]
where
\[(\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)\]
is covariant incompatibility operator in arbitrary dimensions.
[11]:
print("\nRiemann curvature tensor:")
t = 1 / 8
gfGt = GridFunction(gfsigma.space)
errold = None
err = None
with TaskManager():
term_0 = mf.Riemann
variation = -2 * mf.CovInc(sigma) + 0.5 * (
fem.Einsum("ijal,ak->ijkl", mf.Riemann, Inv(gfG) * sigma)
+ fem.Einsum("ijka,al->ijkl", mf.Riemann, Inv(gfG) * sigma)
)
for m in range(8):
t /= 2
gfGt.vec.data = gfG.vec + t * gfsigma.vec
mf_t = dg.RiemannianManifold(gfGt)
term_t = mf_t.Riemann
errold = err
err = sqrt(
Integrate(
InnerProduct(
term_t - term_0 - t * variation, term_t - term_0 - t * variation
),
mesh,
)
)
print(
f"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}"
)
Riemann curvature tensor:
err = 0.1692318844, order = -
err = 0.0478505653, order = 1.822
err = 0.0128254270, order = 1.9
err = 0.0033282435, order = 1.946
err = 0.0008483202, order = 1.972
err = 0.0002141818, order = 1.986
err = 0.0000538127, order = 1.993
err = 0.0000134869, order = 1.996
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
\[\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)\]
[12]:
print("second fundamental form:")
t = 1 / 8
gfGt = GridFunction(gfsigma.space)
errold = None
err = None
Q_eucl = Id(mesh.dim) - OuterProduct(nv, nv)
with TaskManager():
term_0 = Q_eucl * mf.SFF * Q_eucl
variation = (
Q_eucl
* (
0.5 * mf.InnerProduct(sigma, dg.TensorProduct(nv_g, nv_g)) * mf.SFF
+ 0.5
* (
2 * Sym(mf.Contraction(mf.CovDeriv(sigma), nv_g, 1))
- mf.Contraction(mf.CovDeriv(sigma), nv_g, 0)
)
)
* Q_eucl
)
for m in range(8):
t /= 2
gfGt.vec.data = gfG.vec + t * gfsigma.vec
mf_t = dg.RiemannianManifold(gfGt)
term_t = Q_eucl * mf_t.SFF * Q_eucl
errold = err
err = sqrt(
Integrate(
InnerProduct(
term_t - term_0 - t * variation, term_t - term_0 - t * variation
)
* dx(element_boundary=True),
mesh,
)
)
print(
f"err = {err:.10f}, order = {round(log2(errold / err), 3) if errold else '-'}"
)
second fundamental form:
err = 0.8853375262, order = -
err = 0.4661367212, order = 0.925
err = 0.2400109427, order = 0.958
err = 0.1219135777, order = 0.977
err = 0.0614585626, order = 0.988
err = 0.0308580525, order = 0.994
err = 0.0154616510, order = 0.997
err = 0.0077390367, order = 0.998
[ ]: