Distributional Riemann curvature tensor curvature
In this notebook, we discuss the distributional Riemann curvature tensor in arbitrary dimensions. The Riemann curvature tensor stores all the information about the curvature of a Riemannian manifold. From it, we can compute all other curvature quantities, such as the Ricci and Einstein tensor, and the scalar and (in 2D) the Gauss curvature. This also answers the following fundamental question positively: Is it possible to define all intrinsic curvature quantities from a Regge metric? For a detailed derivation and analysis, we refer to Gopalakrishnan, Neunteufel, Schöberl, Wardetzky. Generalizing Riemann curvature to Regge metrics, arXiv (2024)..
Riemann curvature tensor as a nonlinear distribution
The Riemann curvature tensor is a fourth-order tensor with several (skew)-symmetry properties. The distributional (densitized) Riemann curvature tensor acting on specific tensor-valued test functions \(A\in\mathcal{A} is defined as\)$ \begin{align*} \widetilde{\mathfrak{R}\omega}(A) :=\sum_{T\in\mathcal{T}}\int_T\langle \mathfrak{R}|_T,A\rangle\,\omega_T+4\sum_{F\in\mathcal{F}}\int_F\langle[\![I\!I]\!]\,A(\cdot,\nu,\nu,\cdot)\rangle\,\omega_F+4\sum_{E\in\mathcal{E}}\int_E\Theta_E\,A(\mu,\nu,\nu,\mu)\,\omega_E, \end{align*}
\begin{align*} I\!I(X,Y) = g(\nabla_X\nu,Y)=-g(\nabla_XY,\nu),\quad X,Y\in \mathfrak{X}(F), \end{align*}
\begin{align*} A(\cdot,\nu,\nu,\cdot) \text{ is single valued for all $F\in\mathcal{F}$, i.e. } [\![A(\cdot,\nu,\nu,\cdot)]\!]_F=0. \end{align*} $$ From it, we can prove that \(A(\mu,\nu,\nu,\mu)\) is single-valued at all co-dimension 2 boundaries \(E\).
A heuristic motivation of the jump of the second fundamental form as the facet is given as follows: Consider a facet \(F\) that is a subsimplex of two adjacent elements \(T_\pm\) with respective unit normal vectors \(\nu_\pm\) at a point \(p\) in \(F\). Also consider two \(2\)-dimensional submanifolds \(S_\pm\) of \(T_\pm\) passing through \(p\) and transversely intersecting \(F\) along a smooth curve \(\gamma\) in \(F\), such that one tangent direction at \(p\) is parallel to the respective normal \(\nu_{\pm}\). Denote the unit tangent vector along \(\gamma\) at \(p\) by \(X\). Let \(\mathfrak{R}_\pm\) be the respective Riemann curvature tensors of \(T_\pm\). The respective sectional curvatures at \(p\) generated by the orthonormal pair \(X, \nu_\pm\) equal \(\mathfrak{R}_\pm(X, \nu_\pm, \nu_\pm, X)\). It is well-known that these sectional curvatures coincide with the Gauss curvatures of \(S_\pm\). Since we know that the distributional Gauss curvature of the composite \(2\)-manifold \(S_+ \cup S_-\) must have an interface curvature contribution arising from the jump of the geodesic curvatures \(\kappa_\pm \equiv I\!I_\pm(X, X)\) from adjacent elements, we conclude that the difference between \(\mathfrak{R}_\pm(X, \nu_\pm, \nu_\pm, X)\) at \(p\) must equal the difference between \(I\!I_+(X, X)\) and \(I\!I_-(X, X)\). This result can be extended to any \(X\) in the tangent space of \(F\) at \(p\) by varying the 2-manifold. Furthermore, by the symmetry of \(\mathfrak{R}_\pm(X, \nu_\pm, \nu_\pm, Y)\) in \(X\) and \(Y\), the polarization identity implies that the jump of \(\mathfrak{R}_\pm(X, \nu_\pm, \nu_\pm, Y)\) must equal the jump of \(I\!I_\pm(X, Y)\) for any \(X, Y\) in the tangent space of \(F\).
Integral representation and convergence
The distributional Riemann curvature tensor entails an integral representation, which is crucial for numerical analysis. Unfortunately, the test function space \(\mathcal{A}\) also depends on the metric \(g\). We can use the Uhlenbeck trick to overcome this issue, where all quantities get mapped to a metric-independent configuration. Then, the variation is computed, and the result can be mapped back to the current setting. We skip the details and directly present the final integral representation
Taking the variation of the densitized distributional Riemann curvature tensor in the metric in direction \(\sigma\) yields
with
Here, all differential operators are covariant ones with respect to the metric \(g\). Further, \(\mathbb{S}_F(\sigma)=\sigma|_F-\mathrm{tr}(\sigma|_F)g_F\), \((SA)(X,Y,Z,W) =A(X,Z,Y,W)\), and \((L_{\sigma}^{(1)}\mathfrak{R})(X,Y,Z,W) =\mathfrak{R}(L_{\sigma}X,Y,Z,W)\) with \(L_{\sigma}\) the endomorphism defined by \(g(L_{\sigma}X,Y) = \sigma(X,Y)\). Due to the skew-symmetry of \(SA\), the bilinear form \(b_h(g;\sigma,A)\) can be interpreted as a distributional covariant extension of the incompatibility operator in arbitrary dimensions.
As expected, we obtain the same convergence rates as for the distributional scalar curvature in the \(H^{-2}\) norm for \(k\ge 0\) and \(g_h\in \mathrm{Reg}^k\)
Specialization in 3D
In three dimensions, we can use that the curvature operator \(\mathcal{Q}\) is related to the Riemann curvature operator via
where \(X\times Y\) is the covariant cross product defined by \(\omega(X,Y,\cdot)^\sharp\). Further, there holds the following relation between the test functions \(A\in\mathcal{A}\) and the Regge function \(U\in\mathrm{Reg}\)
where \(\hat{\varepsilon}^{ijk}=\sqrt{\det(g)}^{-1}\varepsilon^{ijk}\) is the covariant version of the Levi-Civita symbol. Inserting into the definition of the distributional Riemann curvature tensor and simplifying yields the following distributional curvature operator
where \(\tau_E\) is the \(g\)-unit tangent vector of the edge \(E\) and \((A\times B)^{ij}=\hat{\varepsilon}^{ikl}\hat{\varepsilon}^{jmn}A_{km}B_{ln}\) denotes the tensor cross product. We will use this expression for the numerical example.
Numerical example
In the following example, we compute the distributional Riemann curvature tensor \(\widetilde{\mathfrak{R}\omega}\), more precisely its equivalent curvature operator form \(\widetilde{\mathcal{Q}\omega}\), on a sequence of meshes for a three-dimensional Riemannian manifold and compute the \(H^{-2}\) error.
[1]:
from ngsolve.meshes import MakeStructured3DMesh
from ngsolve import *
import ngsolve.fem as fem
import random as random
from ngsolve.krylovspace import CGSolver
import numpy as np
from ngsolve.webgui import Draw
def mapping(x, y, z):
return (2 * x - 1, 2 * y - 1, 2 * z - 1)
mesh = MakeStructured3DMesh(False, nx=6, ny=6, nz=6, mapping=mapping)
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)
)
G_ex = F.trans * F
def q(x):
return x**2 * (x**2 - 3) ** 2
Qex_xx = 9 * (z**2 - 1) * (y**2 - 1) / (q(x) + q(y) + q(z) + 9)
Qex_xy = 0
Qex_xz = 0
Qex_yy = 9 * (z**2 - 1) * (x**2 - 1) / (q(x) + q(y) + q(z) + 9)
Qex_yz = 0
Qex_zz = 9 * (y**2 - 1) * (x**2 - 1) / (q(x) + q(y) + q(z) + 9)
Q_ex = (
1
/ Det(G_ex)
* CF(
(Qex_xx, Qex_xy, Qex_xz, Qex_xy, Qex_yy, Qex_yz, Qex_xz, Qex_yz, Qex_zz),
dims=(3, 3),
)
)
clipping = {"pnt": (0, 0, 0.01), "z": -1, "function": True}
Draw(Norm(Q_ex), mesh, "Q", clipping=clipping)
Draw(Q_ex[8], mesh, "Qzz", clipping=clipping)
[1]:
<netgen.webgui.WebGLScene at 0x7f9d88097470>
Like for the Einstein tensor, we compute the \(H^{-2}\)-norm of the error \(\|\widetilde{\mathcal{Q}\omega}-\mathcal{Q}\omega\|_{H^{-2}}\) by solving for each component the following fourth-order problem \(\mathrm{div}(\mathrm{div}(\nabla^2 u_{ij}))= (\widetilde{\mathcal{Q}\omega}-\mathcal{Q}\omega)_{ij}\). Then \(\|u\|_{H^2}\approx \|\widetilde{\mathcal{Q}\omega}-\mathcal{Q}\omega\|_{H^{-2}}\). We use the Euclidean Hellan-Herrmann-Johnson (HHJ) method to solve the fourth-order problem as a mixed saddle point problem. Therein, \(\sigma=\nabla^2u_{ij}\) is introduced as an additional unknown. We use hybridization techniques to regain a minimization problem. We use a BDDC preconditioner in combination with static condensation and a CGSolver to solve the problem iteratively. To mitigate approximation errors that spoil the convergence, we use two polynomial degrees more than for \(g_h\) to compute \(u_h\).
[2]:
def ComputeHm2Norm(func_vol, func_bnd, func_bbnd, mesh, order):
fesH = H1(mesh, order=order, dirichlet=".*")
fesS = Discontinuous(HDivDiv(mesh, order=order - 1))
fesHB = NormalFacetFESpace(mesh, order=order - 1, dirichlet=".*")
fes = fesH * fesS * fesHB
(u, sigma, uh), (v, tau, vh) = fes.TnT()
n = specialcf.normal(mesh.dim)
a = BilinearForm(fes, symmetric=True, symmetric_storage=True, condense=True)
a += (
InnerProduct(sigma, tau)
- InnerProduct(tau, u.Operator("hesse"))
- InnerProduct(sigma, v.Operator("hesse"))
) * dx
a += ((Grad(u) - uh) * n * tau[n, n] + (Grad(v) - vh) * n * sigma[n, n]) * dx(
element_boundary=True
)
apre = Preconditioner(a, "bddc")
a.Assemble()
invS = CGSolver(a.mat, apre, printrates="\r", maxiter=600)
if a.condense:
ext = IdentityMatrix() + a.harmonic_extension
inv = a.inner_solve + ext @ invS @ ext.T
else:
inv = invS
# use MultiVector to solve for all right-hand sides at once
rhs = MultiVector(fes.ndof, 6, False)
gf_u_i = GridFunction(fes**6)
for i in range(6):
j = i if i < 3 else (i + 1 if i < 5 else 8)
f = LinearForm(
func_vol[j] * v * dx
+ func_bnd[j] * v * dx(element_boundary=True)
+ func_bbnd[j] * v * dx(element_vb=BBND)
)
f.Assemble()
rhs[i].data = f.vec
gfu_multivec = (inv * rhs).Evaluate()
for i in range(6):
gf_u_i.components[i].vec.data = gfu_multivec[i]
index_symmetric_matrix = [0, 1, 2, 1, 3, 4, 2, 4, 5]
# u
gfu = CF(
tuple([gf_u_i.components[i].components[0] for i in index_symmetric_matrix]),
dims=(3, 3),
)
# Grad(u)
grad_gfu = CF(
tuple(
[Grad(gf_u_i.components[i].components[0]) for i in index_symmetric_matrix]
),
dims=(3, 3, 3),
)
# sigma = hesse u
gf_sigma = CF(
tuple([gf_u_i.components[i].components[1]] for i in index_symmetric_matrix),
dims=(3, 3, 3, 3),
)
err = sqrt(
Integrate(
InnerProduct(gfu, gfu)
+ InnerProduct(grad_gfu, grad_gfu)
+ InnerProduct(gf_sigma, gf_sigma),
mesh,
)
)
return err
The following function computes the \(H^2\)-error of the distributional curvature operator for a given metric and polynomial order.
[3]:
# Euclidean vectors
t = specialcf.tangential(3)
n = -specialcf.normal(3)
bbnd_tang = specialcf.EdgeFaceTangentialVectors(3)
tef1 = bbnd_tang[:, 0]
tef2 = bbnd_tang[:, 1]
# Euclidean tensor-cross product
def TensorCross(A, B):
return fem.Einsum(
"ikl,jmn,km,ln->ij", fem.LeviCivitaSymbol(3), fem.LeviCivitaSymbol(3), A, B
)
# for angle defect
def TetAngle(G):
return acos(G * tef1 * tef2 / sqrt(G * tef1 * tef1) / sqrt(G * tef2 * tef2))
def ComputeCurvature(mesh, G, order):
# Regge metric
fesCC = HCurlCurl(mesh, order=order)
gfG = GridFunction(fesCC)
gfG.Set(G)
# volume forms
vol = sqrt(Det(gfG))
vol_bnd = sqrt(Cof(gfG) * n * n)
vol_bbnd = sqrt(gfG * t * t)
# volume term
func_vol = 1 / vol * gfG.Operator("curvature") - sqrt(Det(G_ex)) * Q_ex
# boundary (facet) term
ng = 1 / sqrt(Inv(gfG) * n * n) * Inv(gfG) * n
func_bnd = (
vol_bnd
/ vol**2
* TensorCross(
OuterProduct(gfG * ng, gfG * ng), gfG.Operator("christoffel") * ng
)
)
# co-dimension 2 (edge) term
tg = 1 / sqrt(gfG * t * t) * t
func_bbnd = vol_bbnd * (TetAngle(Id(3)) - TetAngle(gfG)) * OuterProduct(tg, tg)
err_hm2 = ComputeHm2Norm(func_vol, func_bnd, func_bbnd, mesh, order + 2)
return (err_hm2, fesCC.ndof)
To avoid possible super-convergence due to symmetries in the mesh, we perturb the internal vertices by a random noise.
[4]:
# compute the diameter of the mesh as maximal mesh size h
def CompDiam(mesh):
lengths = np.zeros(mesh.nedge)
for edge in mesh.edges:
pnts = []
for vert in edge.vertices:
pnt = mesh[vert].point
pnts.append((pnt[0], pnt[1], pnt[2]))
lengths[edge.nr] = sqrt(
(pnts[0][0] - pnts[1][0]) ** 2
+ (pnts[0][1] - pnts[1][1]) ** 2
+ (pnts[0][2] - pnts[1][2]) ** 2
)
print(f"max = {max(lengths)}, min = {min(lengths)}, ne = {mesh.ne}")
return max(lengths)
# generate and randomly perturb it, check if the resulting mesh is valid
def GenerateRandomMesh(k, par):
for i in range(20):
mesh = MakeStructured3DMesh(False, nx=2**k, ny=2**k, nz=2**k, mapping=mapping)
maxh = CompDiam(mesh)
for pnts in mesh.ngmesh.Points():
px, py, pz = pnts[0], pnts[1], pnts[2]
if (
abs(px + 1) > 1e-8
and abs(px - 1) > 1e-8
and abs(py + 1) > 1e-8
and abs(py - 1) > 1e-8
and abs(pz + 1) > 1e-8
and abs(pz - 1) > 1e-8
):
px += random.uniform(-maxh / 2**par, maxh / 2**par)
py += random.uniform(-maxh / 2**par, maxh / 2**par)
pz += random.uniform(-maxh / 2**par, maxh / 2**par)
pnts[0] = px
pnts[1] = py
pnts[2] = pz
F = specialcf.JacobianMatrix(3)
# check if the mesh is valid
if Integrate(IfPos(Det(F), 0, -1), mesh) >= 0:
break
else:
raise Exception("Could not generate a valid mesh")
return mesh
Compute the curvature operator on a sequence of meshes. You can play around with the polynomial order. The \(H^{-2}\) error and the number of degrees of freedom are stored.
[5]:
err_hm2 = []
ndof = []
# order of metric approximation >= 0
order = 1
with TaskManager():
for i in range(4):
mesh = GenerateRandomMesh(k=i, par=3)
errm2, dof = ComputeCurvature(mesh, G=G_ex, order=order)
err_hm2.append(errm2)
ndof.append(dof)
max = 3.4641016151377544, min = 2.0, ne = 6
CG converged in 13 iterations to residual 8.100154242394525e-15
CG converged in 10 iterations to residual 5.1811049328140984e-15
CG converged in 10 iterations to residual 2.4789847672814193e-14
CG converged in 13 iterations to residual 7.779964451231337e-15
CG converged in 10 iterations to residual 2.2135948329298652e-14
CG converged in 13 iterations to residual 7.894356995288157e-15
max = 1.7320508075688772, min = 1.0, ne = 48
CG converged in 39 iterations to residual 1.6892514820107565e-14
CG converged in 39 iterations to residual 7.30643286145417e-15
CG converged in 39 iterations to residual 7.466115196674338e-15
CG converged in 39 iterations to residual 1.58547943020992e-14
CG converged in 40 iterations to residual 4.499841742458574e-15
CG converged in 40 iterations to residual 9.226385071067919e-15
max = 0.8660254037844386, min = 0.5, ne = 384
CG converged in 61 iterations to residual 9.535287545798354e-15
CG converged in 61 iterations to residual 3.4821936217553516e-15
CG converged in 61 iterations to residual 3.4713369819476537e-15
CG converged in 61 iterations to residual 9.245779857421073e-15
CG converged in 61 iterations to residual 3.7466545734287266e-15
CG converged in 61 iterations to residual 8.432437168106885e-15
max = 0.4330127018922193, min = 0.25, ne = 3072
CG converged in 94 iterations to residual 3.0227376995328676e-15
CG converged in 94 iterations to residual 1.3318871911382708e-15
CG converged in 94 iterations to residual 1.3649060233369468e-15
CG converged in 94 iterations to residual 3.062590972604862e-15
CG converged in 94 iterations to residual 1.3085354596455395e-15
CG converged in 94 iterations to residual 3.0904846192067432e-15
Plot the results together with reference lines for linear, quadratic, and cubic convergence order.
[6]:
import matplotlib.pyplot as plt
plt.plot(
ndof,
err_hm2,
"-x",
label=r"$\|\widetilde{\mathcal{Q}\omega}-\mathcal{Q}_{\mathrm{ex}}\omega_{\mathrm{ex}}\|_{H^{-2}}$",
)
plt.plot(ndof, [2 * dof ** (-1 / 3) for dof in ndof], "-", color="k", label="$O(h)$")
plt.plot(ndof, [5 * dof ** (-2 / 3) for dof in ndof], "--", color="k", label="$O(h^2)$")
plt.plot(
ndof, [50 * dof ** (-3 / 3) for dof in ndof], "-.", color="k", label="$O(h^3)$"
)
plt.yscale("log")
plt.xscale("log")
plt.legend()
plt.title(f"order = {order}")
plt.show()
[ ]: