Curvature quantities on Riemannian manifolds

In this notebook, we discuss the definition and computation of different curvature quantities of Riemannian manifolds and their relations, see this notebook for a brief introduction to Riemannian manifolds. We focus on the Riemann curvature tensor, curvature operator, Ricci curvature tensor, scalar curvature, and (in 3D) the Einstein tensor.

We use two- and three-dimensional numerical examples to verify the identities. We consider the Poincaré disk in 2D and the 3-sphere in 3D. For the numerical computations, we interpolate the exact metric into the Regge finite elements.

[1]:

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import *
import ngsdiffgeo as dg
from ngsolve.fem import Einsum, LeviCivitaSymbol
from time import time

mesh2d = Mesh(
    OCCGeometry(Circle((0, 0), 0.7).Face(), dim=2).GenerateMesh(maxh=0.05)
).Curve(3)
pcd = dg.PoincareDisk()
Draw(Norm(pcd.metric), mesh2d)

eps = 2e-1
mesh3d = Mesh(
    OCCGeometry(Box((eps, eps, eps), (pi - eps, pi - eps, 2 * pi - eps))).GenerateMesh(
        maxh=0.4
    )
)
sphere3 = dg.Sphere3()
Draw(Norm(sphere3.metric), mesh3d, order=3)


# Euclidean Levi-Civita permuting symbols in 2D and 3D
Eps2d = LeviCivitaSymbol(2)
Eps3d = LeviCivitaSymbol(3)

# Interpolate exact metric to Regge space
g2d = GridFunction(HCurlCurl(mesh2d, order=4))
g3d = GridFunction(HCurlCurl(mesh3d, order=3))

with TaskManager():
    g2d.Set(pcd.metric)
    g3d.Set(sphere3.metric)

# Define Christoffel symbols of first and second kind
chr1_2d = g2d.Operator("christoffel")
chr2_2d = g2d.Operator("christoffel2")
chr1_3d = g3d.Operator("christoffel")
chr2_3d = g3d.Operator("christoffel2")

mf2d = dg.RiemannianManifold(g2d)
mf3d = dg.RiemannianManifold(g3d)

<netgen.webgui.WebGuiWidget>: Run `pip install anywidget` to enable WebGuiWidget in Jupyter

<netgen.webgui.WebGuiWidget>: Run `pip install anywidget` to enable WebGuiWidget in Jupyter

Riemann curvature tensor

We start with the fourth-order Riemann curvature tensor \(\mathfrak{R}\). Its coordinate expression is

\[\begin{align*} \mathfrak{R}_{ijkl} = \partial_i \Gamma_{jkl}-\partial_j \Gamma_{ikl} + \Gamma_{ik}^p\Gamma_{jlp}- \Gamma_{jk}^p\Gamma_{ilp}. \end{align*}\]

First, we note that the first two terms can be rewritten using the incompatibility operator \(\mathrm{inc}(g)= \widehat\varepsilon^{ij}\widehat\varepsilon^{kl}\partial_i\partial_k g_{jl}\) in 2D and \(\mathrm{inc}(g)^{ij}= \widehat\varepsilon^{ikl}\widehat\varepsilon^{jmn}\partial_k\partial_m g_{ln}\) in 3D. Here \(\widehat\varepsilon\) denotes the alternating permutation symbol, as used by LeviCivitaSymbol in the code. There holds

\[\begin{split}\begin{align*} \partial_i \Gamma_{jkl}-\partial_j \Gamma_{ikl} = \begin{cases}\frac{1}{2}\widehat\varepsilon_{ij}\widehat\varepsilon_{kl}\mathrm{inc}(g), & \text{in 2D},\\ \frac{1}{2}\widehat\varepsilon_{ija}\widehat\varepsilon_{klb}\mathrm{inc}(g)^{ab}, & \text{in 3D}. \end{cases} \end{align*}\end{split}\]

We compare the Riemann differential operator from the Regge finite element with its coordinate expression for the same interpolated metric \(g_h\), and separately with the exact curvature of the smooth metric. The two expressions for \(g_h\) agree up to numerical roundoff. Their error relative to the exact curvature includes metric-interpolation error and can be substantial on the coarse 3D mesh. The later relative-error values should not be read as roundoff errors. Mesh refinement reduces this approximation error. We will learn later how to approximate the curvature more accurately from the metric tensor.

[2]:

# NGSolve Riemann curvature tensor
Riem1_2d = g2d.Operator("Riemann")
# manually compute Riemann curvature tensor
Riem2_2d = (
    0.5 * g2d.Operator("inc") * Einsum("ij,kl->ijkl", Eps2d, Eps2d)
    + Einsum("ikp,jlp->ijkl", chr2_2d, chr1_2d)
    - Einsum("jkp,ilp->ijkl", chr2_2d, chr1_2d)
)
# exact Riemann curvature tensor
Riem_ex_2d = pcd.Riemann

# The same in 3D
Riem1_3d = g3d.Operator("Riemann")
Riem2_3d = (
    0.5 * Einsum("ija,klb,ab->ijkl", Eps3d, Eps3d, g3d.Operator("inc"))
    + Einsum("ikp,jlp->ijkl", chr2_3d, chr1_3d)
    - Einsum("jkp,ilp->ijkl", chr2_3d, chr1_3d)
)
Riem_ex_3d = sphere3.Riemann

# check that the expressions are equal
with TaskManager():
    print(
        f"err Riemann curvature tensor 2d = {sqrt(Integrate(InnerProduct(Riem1_2d - Riem2_2d, Riem1_2d - Riem2_2d), mesh2d)):.8f}"
    )
    print(
        f"rel err Riemann curvature tensor 2d exact = {sqrt(Integrate(InnerProduct(Riem1_2d - Riem_ex_2d, Riem1_2d - Riem_ex_2d), mesh2d)) / sqrt(Integrate(InnerProduct(Riem_ex_2d, Riem_ex_2d), mesh2d)):.8f}"
    )
    print(
        f"err Riemann curvature tensor 3d = {sqrt(Integrate(InnerProduct(Riem1_3d - Riem2_3d, Riem1_3d - Riem2_3d), mesh3d)):.8f}"
    )

    print(
        f"rel err Riemann curvature tensor 3d exact = {sqrt(Integrate(InnerProduct(Riem1_3d - Riem_ex_3d, Riem1_3d - Riem_ex_3d), mesh3d)) / sqrt(Integrate(InnerProduct(Riem_ex_3d, Riem_ex_3d), mesh3d)):.8f}"
    )

err Riemann curvature tensor 2d = 0.00000000
rel err Riemann curvature tensor 2d exact = 0.00188534
err Riemann curvature tensor 3d = 0.00000000
rel err Riemann curvature tensor 3d exact = 0.03248129

Curvature operator

Next, we look at the curvature operator \(\mathcal{Q}\), which incorporates the skew-symmetry property of the Riemann curvature tensor. It is defined via the identity

\[\begin{align*} \mathcal{Q}(X^\flat\wedge Y^\flat, W^\flat \wedge Z^\flat) = \mathfrak{R}(X,Y,Z,W),\qquad \forall X,Y,Z,W\in\mathfrak{X}(M). \end{align*}\]

Note the convention that the \(Z\) and \(W\) arguments are swapped and that \(X^\flat\wedge Y^\flat \in \Lambda^2(M)\) is a 2-form. In two dimensions the curvature operator can be identified with a scalar function as \(\mathcal{Q}\,\omega(X,Y)\,\omega(W,Z) = \mathcal{Q}(X^\flat\wedge Y^\flat, W^\flat\wedge Z^\flat)\) because 2-forms and the volume form are linearly dependent. Thus, in coordinates, there holds

\[\begin{align*} \mathcal{Q} = -\frac{1}{4}\varepsilon^{kl}\varepsilon^{mn}\mathfrak{R}_{klmn} = -\frac{1}{4\,\det g}\widehat\varepsilon^{kl}\widehat\varepsilon^{mn}\mathfrak{R}_{klmn}, \end{align*}\]

where \(\varepsilon^{kl}= \frac{1}{\sqrt{\det g}}\widehat\varepsilon^{kl}\) are the contravariant components of the Levi-Civita tensor. In three dimensions, \(X^\flat\wedge Y^\flat\) can be identified with the cross product \((X\times Y)^\flat\). The curvature operator becomes a symmetric operator acting on two 1-forms and reads in coordinates

\[\begin{align*} \mathcal{Q}^{ij} = -\frac{1}{4}\varepsilon^{ikl}\varepsilon^{jmn}\mathfrak{R}_{klmn} = -\frac{1}{4\,\det g}\widehat\varepsilon^{ikl}\widehat\varepsilon^{jmn}\mathfrak{R}_{klmn}. \end{align*}\]

Using the coordinate expression of the Riemann curvature tensor and the alternating-symbol contraction identity \(\widehat\varepsilon^{ija}\widehat\varepsilon_{kla}=\delta^i_k\delta^j_l-\delta^i_l\delta^j_k\), we obtain in 2D and 3D

\[\begin{align*} \mathcal{Q}=-\frac{1}{2}\varepsilon^{kl}\varepsilon^{mn}\left(\partial_k\partial_m g_{ln}+\Gamma^q_{km}\Gamma_{lnq}\right),\qquad \mathcal{Q}^{ij}=-\frac{1}{2}\varepsilon^{ikl}\varepsilon^{jmn}\left(\partial_k\partial_m g_{ln}+\Gamma^q_{km}\Gamma_{lnq}\right). \end{align*}\]

Conversely, we can express the Riemann curvature tensor in terms of the curvature operator

\[\begin{split}\begin{align*} \mathfrak{R}_{ijkl}=\begin{cases} -\varepsilon_{ij}\varepsilon_{kl}\mathcal{Q}, &\text{in 2D}\\ -\varepsilon_{ija}\varepsilon_{klb}\mathcal{Q}^{ab}&\text{in 3D}.\end{cases} \end{align*}\end{split}\]

Here \(\varepsilon_{ij}=\sqrt{\det g}\,\widehat\varepsilon_{ij}\) in 2D, with the analogous relation in 3D. The minus signs follow from the argument order \(\mathfrak{R}(X,Y,Z,W)=\mathcal{Q}(X^\flat\wedge Y^\flat,W^\flat\wedge Z^\flat)\).

WARNING: The curvature differential operator .Operator(“curvature”) of a Regge object returns \(\det g\,\mathcal{Q}\) instead of \(\mathcal{Q}\)! Thus, we need to divide by the determinant.

mf.Curvature returns the normalized geometric operator for both ordinary coefficient-function and Regge metrics. In 2D it equals mf.Gauss and half the scalar curvature. In 3D it is contravariant and satisfies \(\mathrm{Einstein}=-gQg\). NGSolve’s native g.Operator("curvature") is unnormalized and must still be divided by Det(g).

Compatibility note: earlier NGSDiffGeo versions returned the unnormalized operator for Regge inputs through mf.Curvature. Remove any additional division by Det(g) when migrating such calls. Direct native-operator expressions keep their existing normalization.

[3]:

Q1_2d = mf2d.Curvature
Q2_2d = (
    -1 / (4 * Det(g2d)) * Einsum("ij,kl,ijkl", Eps2d, Eps2d, g2d.Operator("Riemann"))
)
Q3_2d = (
    -0.5
    / Det(g2d)
    * (g2d.Operator("inc") + Einsum("kl,mn,kmq,lnq", Eps2d, Eps2d, chr2_2d, chr1_2d))
)
Q_ex_2d = pcd.curvature

Q1_3d = mf3d.Curvature
Q2_3d = (
    -1
    / (4 * Det(g3d))
    * Einsum("ikl,jmn,klmn->ij", Eps3d, Eps3d, g3d.Operator("Riemann"))
)
Q3_3d = (
    -0.5
    / Det(g3d)
    * (
        g3d.Operator("inc")
        + Einsum("ikl,jmn,kmq,lnq->ij", Eps3d, Eps3d, chr2_3d, chr1_3d)
    )
)
Q_ex_3d = sphere3.curvature

with TaskManager():
    print(
        f"err curvature operator 2d 1-2 = {sqrt(Integrate((Q1_2d - Q2_2d) ** 2, mesh2d)):.8f}"
    )
    print(
        f"err curvature operator 2d 1-3 = {sqrt(Integrate((Q1_2d - Q3_2d) ** 2, mesh2d)):.8f}"
    )
    print(
        f"rel err curvature operator 2d exact = {sqrt(Integrate((Q1_2d - Q_ex_2d) ** 2, mesh2d)) / sqrt(Integrate(Q_ex_2d**2, mesh2d)):.8f}"
    )

    print(
        f"err curvature operator 3d 1-2 = {sqrt(Integrate(InnerProduct(Q1_3d - Q2_3d, Q1_3d - Q2_3d), mesh3d)):.8f}"
    )
    print(
        f"err curvature operator 3d 1-3 = {sqrt(Integrate(InnerProduct(Q1_3d - Q3_3d, Q1_3d - Q3_3d), mesh3d)):.8f}"
    )
    print(
        f"rel err curvature operator 3d exact = {sqrt(Integrate(InnerProduct(Q1_3d - Q_ex_3d, Q1_3d - Q_ex_3d), mesh3d)) / sqrt(Integrate(InnerProduct(Q_ex_3d, Q_ex_3d), mesh3d)):.8f}"
    )

err curvature operator 2d 1-2 = 0.00000000
err curvature operator 2d 1-3 = 0.00000000
rel err curvature operator 2d exact = 0.00100136
err curvature operator 3d 1-2 = 0.00000000
err curvature operator 3d 1-3 = 0.00000000
rel err curvature operator 3d exact = 0.46355484

Ricci curvature tensor

The Ricci curvature tensor is defined by contracting two of the indices of the Riemann curvature tensor (i.e., taking a trace)

\[\begin{align*} \mathrm{Ric}_{ij} = \mathfrak{R}_{kijl}g^{kl}, \end{align*}\]

In terms of the curvature operator, we have

\[\begin{split}\begin{align*} \mathrm{Ric} = \begin{cases} \mathcal{Q}\,g, &\text{in 2D}\\ g^{-1}\times \mathcal{Q}, &\text{in 3D}\end{cases}, \end{align*}\end{split}\]

where \((A \times B)_{ij}= \varepsilon_{ikl}\varepsilon_{jmn}A^{km}B^{ln}\) is the tensor cross product.

[4]:

Ric1_2d = g2d.Operator("Ricci")
Ric2_2d = 1 / Det(g2d) * g2d.Operator("curvature") * g2d
Ric_ex_2d = pcd.Ricci

Ric1_3d = g3d.Operator("Ricci")
Ric2_3d = Det(g3d) * Einsum(
    "ikl,jmn,km,ln->ij",
    Eps3d,
    Eps3d,
    Inv(g3d),
    1 / Det(g3d) * g3d.Operator("curvature"),
)
Ric_ex_3d = sphere3.Ricci

with TaskManager():
    print(
        f"err Ricci curvature tensor 2d = {sqrt(Integrate(InnerProduct(Ric1_2d - Ric2_2d, Ric1_2d - Ric2_2d), mesh2d)):.8f}"
    )
    print(
        f"rel err Ricci curvature tensor 2d exact = {sqrt(Integrate(InnerProduct(Ric1_2d - Ric_ex_2d, Ric1_2d - Ric_ex_2d), mesh2d)) / sqrt(Integrate(InnerProduct(Ric_ex_2d, Ric_ex_2d), mesh2d)):.8f}"
    )
    print(
        f"err Ricci curvature tensor 3d = {sqrt(Integrate(InnerProduct(Ric1_3d - Ric2_3d, Ric1_3d - Ric2_3d), mesh3d)):.8f}"
    )
    print(
        f"rel err Ricci curvature tensor 3d exact = {sqrt(Integrate(InnerProduct(Ric1_3d - Ric_ex_3d, Ric1_3d - Ric_ex_3d), mesh3d)) / sqrt(Integrate(InnerProduct(Ric_ex_3d, Ric_ex_3d), mesh3d)):.8f}"
    )

err Ricci curvature tensor 2d = 0.00000000
rel err Ricci curvature tensor 2d exact = 0.00151908
err Ricci curvature tensor 3d = 0.00000000
rel err Ricci curvature tensor 3d exact = 0.21799012

Scalar curvature

The scalar curvature is given by contracting the Riemann curvature tensor twice (or equivalently taking the trace of the Ricci curvature tensor)

\[\begin{align*} S = \mathrm{Ric}_{ij} g^{ij} = \mathfrak{R}_{kijl}g^{kl}g^{ij}. \end{align*}\]

In two dimensions, the Gauss curvature is given for linear independent \(X,Y\) by

\[\begin{align*} K = \frac{\mathcal{R}(X,Y,Y,X)}{g(X,X)g(Y,Y)-g(X,Y)^2}, \end{align*}\]

which is independent of the choice of \(X,Y\). Further, it is related to the scalar curvature by \(S=2K\).

In two dimensions, the Gauss curvature \(K\) coincides with the curvature operator \(\mathcal{Q}\), \(K=\mathcal{Q}\). In 3D, the scalar curvature can be computed from the curvature operator by taking its trace

\[\begin{align*} S = 2 \mathrm{tr}_g(\mathcal{Q}) = 2 g_{ij}\mathcal{Q}^{ij}. \end{align*}\]

[5]:

S1_2d = g2d.Operator("scalar")
S2_2d = Einsum("il,jk,ijkl", Inv(g2d), Inv(g2d), g2d.Operator("Riemann"))
S3_2d = 2 / Det(g2d) * g2d.Operator("curvature")
S_ex_2d = pcd.scalar

S1_3d = g3d.Operator("scalar")
S2_3d = Einsum("il,jk,ijkl", Inv(g3d), Inv(g3d), g3d.Operator("Riemann"))
S3_3d = 2 / Det(g3d) * InnerProduct(g3d, g3d.Operator("curvature"))
S_ex_3d = sphere3.scalar

with TaskManager():
    print(
        f"err scalar curvature 2d 1-2 = {sqrt(Integrate((S1_2d - S2_2d) ** 2, mesh2d)):.8f}"
    )
    print(
        f"err scalar curvature 2d 1-3 = {sqrt(Integrate((S1_2d - S3_2d) ** 2, mesh2d)):.8f}"
    )
    print(
        f"rel err scalar curvature 2d exact = {sqrt(Integrate((S1_2d - S_ex_2d) ** 2, mesh2d)) / sqrt(Integrate(S_ex_2d**2, mesh2d)):.8f}"
    )

    print(
        f"err scalar curvature 3d 1-2 = {sqrt(Integrate((S1_3d - S2_3d) ** 2, mesh3d)):.8f}"
    )
    print(
        f"err scalar curvature 3d 1-3 = {sqrt(Integrate((S1_3d - S3_3d) ** 2, mesh3d)):.8f}"
    )
    print(
        f"rel err scalar curvature 3d exact = {sqrt(Integrate((S1_3d - S_ex_3d) ** 2, mesh3d)) / sqrt(Integrate(S_ex_3d**2, mesh3d)):.8f}"
    )

err scalar curvature 2d 1-2 = 0.00000000
err scalar curvature 2d 1-3 = 0.00000000
rel err scalar curvature 2d exact = 0.00100136
err scalar curvature 3d 1-2 = 0.00000000
err scalar curvature 3d 1-3 = 0.00000000
rel err scalar curvature 3d exact = 0.76544376

Einstein tensor

The Einstein tensor is defined by

\[\begin{align*} G_{ij}=\mathrm{Ric}_{ij}-\frac{1}{2}S\,g_{ij},\qquad G = J\mathrm{Ric}, \end{align*}\]

where \(J\sigma = \sigma-\frac{1}{2}\mathrm{tr}(\sigma)\,g\). In three dimensions, it is the negative of the double lowering of the curvature operator

\[\begin{align*} G = -\mathcal{Q}^{\flat\flat},\qquad G_{ij} = -g_{ik}\mathcal{Q}^{kl}g_{lj}. \end{align*}\]

[6]:

G1_3d = g3d.Operator("Einstein")
G2_3d = g3d.Operator("Ricci") - 0.5 * g3d.Operator("scalar") * g3d
G3_3d = -g3d * (1 / Det(g3d) * g3d.Operator("curvature")) * g3d
G_ex_3d = sphere3.Einstein

with TaskManager():
    print(
        f"err Einstein tensor 3d 1-2 = {sqrt(Integrate(InnerProduct(G1_3d - G2_3d, G1_3d - G2_3d), mesh3d)):.8f}"
    )
    print(
        f"err Einstein tensor 3d 1-3 = {sqrt(Integrate(InnerProduct(G1_3d - G3_3d, G1_3d - G3_3d), mesh3d)):.8f}"
    )
    print(
        f"rel err Einstein tensor 3d exact = {sqrt(Integrate(InnerProduct(G1_3d - G_ex_3d, G1_3d - G_ex_3d), mesh3d)) / sqrt(Integrate(InnerProduct(G_ex_3d, G_ex_3d), mesh3d)):.8f}"
    )

err Einstein tensor 3d 1-2 = 0.00000000
err Einstein tensor 3d 1-3 = 0.00000000
rel err Einstein tensor 3d exact = 1.80833919

[ ]: