Covariant Derivatives

In this notebook, we show how to differentiate scalar, vector, 1-form, and arbitrary tensor fields using the unique Levi-Civita connection of a Riemannian manifold. A brief introduction to Riemannian manifolds is given here. The resulting differential operators are denoted as covariant derivatives, such as the covariant divergence and covariant curl. They simplify to the usual derivatives in the Euclidean setting. In the presence of curvature, the coordinate expressions involve additional terms using the Christoffel symbols of the Levi-Civita connection.

First, we define several fields and discuss their inner products. Then, the covariant derivatives are defined and explained in coordinates. We will see that integration by parts formulae also hold in the general Riemannian setting, as do (some) differential identities known from the Euclidean case.

[1]:

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

TOL = 1e-7

Let \((M,g)\) be a Riemannian manifold with metric tensor \(g\) on a 2D manifold \(M\). Charts induce local coordinates \(x^i\) in the local basis \(\partial_i\) with their corresponding dual 1-forms \(dx^i\). We define the volume forms in the local coordinates on the elements and edges, \begin{align*} \omega = \sqrt{\det(g_{ij})}\,dx^1\wedge\dots\wedge dx^N,\qquad \omega_{\partial M} = \star n^\flat, \end{align*} where \(n^\flat\) is the 1-form associated with the \(g\)-normalized normal vector \(n\) and \(\star\) the Hodge-star operator.

[2]:

metric = dg.CigarSoliton().metric
christoffel1 = dg.CigarSoliton().chr1
christoffel2 = dg.CigarSoliton().chr2
mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))
mf = dg.RiemannianManifold(metric)
omega_T = mf.VolumeForm(VOL)
omega_S = mf.VolumeForm(BND)

Scalar, vector, 1-form, and tensor fields

We start by defining some scalar, vector, 1-form, and tensor fields

\begin{align*} &f,g,h\in C^\infty(M),\qquad X = X^i\,\partial_i\in\mathfrak{X}(M),\qquad \alpha = \alpha_i\,dx^i\in \Lambda^1(M),\\ & A= A^{ij}\partial_i\otimes \partial_j\in\mathcal{T}^0_2(M),\quad B = \alpha\otimes X = \alpha_iX^j\,dx^i\otimes \partial_j\in\mathcal{T}^1_1(M). \end{align*} We use a string to indicate whether the indices of a general tensor field are covariant or contravariant, where 1 denotes covariant indices and 0 denotes contravariant indices.

[3]:

f = dg.ScalarField(CF(x**2 * y - 0.1 * y * x), dim=2)
g = dg.ScalarField(CF(x * y**2 + 0.1 * y * x - 0.73 * x), dim=2)
h = dg.TensorProduct(f, g)

Z = dg.OneForm(dg.GradCF(f, 2))

X = dg.VectorField(CF((10 * x * y**3 - x**2, y**4 * x - y)))
Y = dg.VectorField(CF((y**2 * x - 0.1 * y * x**2, 10 * x * y**3 + x**2 - y)))

alpha = dg.OneForm(CF((sin(x * y), cos(x) * y**2)))
beta = dg.OneForm(
    CF((0.2 * y * x**2 + 0.37 * x**3, 2 * x**2 * y**2 + 0.1 * x**2 - 1.34 * y))
)

A = dg.TensorField(
    CF((10 * x * y**3 - x**2, y**4 * x - y, sin(x * y), cos(x) * y**2), dims=(2, 2)),
    "00",
)
B = dg.TensorProduct(alpha, X)

C = dg.TensorProduct(alpha, beta)

Inner product

The \(g\)-inner product of scalar \(f\), \(g\), vector \(X\), \(Y\), and 1-form \(\alpha\), \(\beta\) fields is given by

\begin{align*} g(f,g) = f\,g,\qquad g(X,Y)= g_{ij}X^iY^j,\qquad g(\alpha,\beta) = g^{ij}\alpha_i\beta_j, \qquad g(\alpha,X) = \alpha_i\beta^i, \end{align*} where we used the Einstein summation convention of repeated indices, e.g.

\begin{align*} g_{ij}X^iY^j:=\sum_{i,j}g_{ij}X^iY^j. \end{align*}

For two tensor fields of the form \(A=X_1 \otimes \dots \otimes X_k \otimes \alpha_1\otimes\dots\otimes \alpha_l\) and \(B=Y_1\otimes\dots\otimes Y_m\otimes \beta_1\otimes\dots\otimes \beta_n\) of the same rank \(r=k+l=m+n\) the definition of the inner product of vector and 1-form fields directly extends to \begin{align*} g(A,B) := g(X_1,Y_1)\cdots g(\alpha_l,\beta_n). \end{align*}

The InnerProduct method of the Riemannian manifold class recognizes the type of tensor fields and employs the correct inner product.

[4]:

left = Integrate(mf.InnerProduct(f, g), mesh)
right = Integrate(f * g, mesh)
print(f"{left:.6f} = {right:.6f}")
assert abs(left - right) < TOL

left = Integrate(mf.InnerProduct(X, Y), mesh)
right = Integrate(metric[X, Y], mesh)
print(f"{left:.6f} = {right:.6f}")
assert abs(left - right) < TOL

left = Integrate(mf.InnerProduct(alpha, Y), mesh)
right = Integrate(InnerProduct(alpha, Y), mesh)
print(f"{left:.6f} = {right:.6f}")
assert abs(left - right) < TOL

left = Integrate(mf.InnerProduct(alpha, beta), mesh)
right = Integrate(Inv(metric)[alpha, beta], mesh)
print(f"{left:.6f} = {right:.6f}")
assert abs(left - right) < TOL

left = Integrate(mf.InnerProduct(A, A), mesh)
right = Integrate(InnerProduct(metric * A * metric, A), mesh)
print(f"{left:.6f} = {right:.6f}")
assert abs(left - right) < TOL

left = Integrate(mf.InnerProduct(A, B), mesh)
right = Integrate(InnerProduct(metric * A, B), mesh)
print(f"{left:.6f} = {right:.6f}")
assert abs(left - right) < TOL

-0.017694 = -0.017694
-0.019774 = -0.019774
0.575764 = 0.575764
-0.203729 = -0.203729
0.838426 = 0.838426
1.010107 = 1.010107

Covariant differential operators

Covariant derivative

Covariant derivatives enable the differentiation of tensor fields on Riemannian manifolds. It takes a \((k,l)\)-tensor field \(\mathcal{T}^k_l(M)\) acting on \(k\) vector fields and \(l\) covectors (1-forms) and returns a \((k+1,l)\)-tensor \(\mathcal{T}^{k+1}l(M)\), i.e. for \(A\in \mathcal{T}^k_l(M)\) and \(X,Y_1,\dots, Y_k\in \mathfrak{X}(M)\), \(\alpha_1,\dots, \alpha_l\) yields \(\nabla A\in \mathcal{T}^{k+1}l(M)\) with

\begin{align*} &(\nabla A)(X,Y_1,\dots, Y_k,\alpha_1,\dots, \alpha_l)= (\nabla_XA)(Y_1,\dots. Y_k,\alpha_1,\dots, \alpha_l)\\ & = X(A(Y_1,\dots, Y_k,\alpha_1,\dots, \alpha_l)) - A(\nabla_XY_1,\dots Y_k,\alpha_1,\dots \alpha_l) -\cdots - A(Y_1,\dots, Y_k,\alpha_1,\dots, \nabla_X\alpha_l), \end{align*}

where the last equality is the Leibnitz rule, which can be used to define the covariant derivative for arbitrary tensors from the differentiation of vector and 1-form fields.

In the flat case, these operators simplify to their Euclidean counterparts. Therefore, in coordinates, they are expressed by standard derivatives with additional terms involving the metric and its derivatives. Particularly the Christoffel symbols of first \(\Gamma_{ijk}\) and second \(\Gamma_{ij}^k\) kind,

\begin{align*} \Gamma_{ijk}=\frac{1}{2}\left(\partial_i g_{jk}+\partial_jg_{ik}-\partial_kg_{ij}\right),\qquad \Gamma_{ij}^k = g^{kl}\Gamma_{ijl} =\frac{1}{2}g^{kl}\left(\partial_i g_{jl}+\partial_jg_{il}-\partial_lg_{ij}\right). \end{align*}

Further, the covariant derivative is compatible with the metric \(g\), and thus its inverse and determinant (volume forms), by construction (the unique Levi-Civita connection). That means

\begin{align*} \nabla g = 0,\quad \nabla g^{-1} = 0,\quad \nabla \omega=0. \end{align*} However, expressing the volume form in the local basis \(\omega=\sqrt{\det (g_{ij})}\,dx^1\wedge\dots\wedge dx^N\) the covariant derivative of the determinant is in general not zero \begin{align*} \nabla \sqrt{\det g}\neq 0. \end{align*}

For a scalar field \(s\in C^\infty(M)\), the covariant derivative coincides with the exterior derivative (and Lie derivative), yielding the 1-form \(\nabla s= ds\in \Lambda^1(M)\). The covariant gradient is obtained by the musical \(\sharp\) operator giving a vector field \begin{align*} \nabla s = d s = \partial_i s\,dx^i,\qquad \mathrm{grad}(s) = (\nabla s)^\sharp = g^{ij}\partial_j s\,\partial_i. \end{align*} Thus, \(\nabla_X s = (ds)(X) = X(s) = X^i\partial_is\) for a vector field \(X= X^i\partial_i\in \mathfrak{X}(M)\).

For vector field \(v\in \mathfrak{X}(M)\) and a 1-form \(\alpha\in \Lambda^1(M)\) the covariant derivative is given in coordinates by \begin{align*} &\nabla v = (\partial_i v^j + \Gamma^j_{ik}v^k)dx^i\otimes \partial_j,\qquad \nabla\alpha = (\partial_i \alpha_j - \Gamma_{ij}^k\alpha_k)dx^i\otimes dx^j. \end{align*} Note the different signs in front of the Christoffel symbols, whether they are contracted with a vector field or 1-form. Using the Leibnitz rule above, the covariant derivative of an arbitrary order tensor \(A\in \mathcal{T}^k_l(M)\) is now easily given \begin{align*} &\nabla A = \left(\partial_i A_{r_1\dots r_k}^{s_1\dots s_l} -\sum_{j=1}^k\Gamma^{t}_{ir_j}A^{s_1\dots s_l}_{r_1\dots t\dots r_k}+ \sum_{j=1}^l\Gamma^{s_j}_{it}A^{s_1\dots t\dots s_l}_{r_1\dots r_k}\right)dx^i\otimes dx^{r_1}\otimes\dots \otimes dx^{r_k}\otimes \partial_{s_1}\otimes\dots\otimes \partial_{s_l}. \end{align*}

The covariant derivative method of the Riemannian manifold class takes the tensor’s signature into account. Further, a so-called GradCF (gradient coefficient function) is used to compute the derivative of an arbitrary tensor numerically. It differentiates classical CoefficientFunctions and GridFunctions. Differentiation of trial and test functions is WIP.

[5]:

# scalar field
left = Integrate(mf.CovDeriv(f), mesh).Norm()
right = Integrate(CF((f.Diff(x), f.Diff(y))), mesh).Norm()
print(f"{left:.6f} = {right:.6f}")
assert abs(left - right) < TOL

# vector field
cov_grad_X = CF((X.Diff(x), X.Diff(y)), dims=(2, 2)) + Einsum(
    "ikj,k->ij", christoffel2, X
)
left = Integrate(mf.CovDeriv(X), mesh).Norm()
right = Integrate(cov_grad_X, mesh).Norm()
print(f"{left:.6f} = {right:.6f}")
assert abs(left - right) < TOL

# 1-form
cov_grad_alpha = CF((alpha.Diff(x), alpha.Diff(y)), dims=(2, 2)) - Einsum(
    "ijk,k->ij", christoffel2, alpha
)
left = Integrate(mf.CovDeriv(alpha), mesh).Norm()
right = Integrate(cov_grad_alpha, mesh).Norm()
print(f"{left:.6f} = {right:.6f}")
assert abs(left - right) < TOL

# second order tensor field
cov_grad_B = (
    CF((B.Diff(x), B.Diff(y)), dims=(2, 2, 2))
    - Einsum("ijl,lk->ijk", christoffel2, B)
    + Einsum("ilk,jl->ijk", christoffel2, B)
)
left = Integrate(mf.CovDeriv(B), mesh).Norm()
right = Integrate(cov_grad_B, mesh).Norm()
print(f"{left:.6f} = {right:.6f}")
assert abs(left - right) < TOL

# covariant derivative of metric is zero (Levi-Civita connection is metric-compatible)
should_be_zero = Integrate(mf.CovDeriv(mf.G), mesh).Norm()
print(f"{should_be_zero:.6f} = 0")
assert should_be_zero < TOL
should_be_zero = Integrate(mf.CovDeriv(mf.G_inv), mesh).Norm()
print(f"{should_be_zero:.6f} = 0")
assert should_be_zero < TOL


# covariant derivative of volume form expressed as determinant in local coordinates
# is in general not zero
print(f"{Integrate(mf.CovDeriv(dg.ScalarField(omega_T, dim=2)), mesh).Norm():.6f} != 0")

0.531769 = 0.531769
4.814188 = 4.814188
1.144794 = 1.144794
4.772021 = 4.772021
0.000000 = 0
0.000000 = 0
0.495241 != 0

Covariant divergence

The covariant divergence of an arbitrary tensor (of order at least 1) is defined as first taking the covariant derivative and then contracting the first two indices, i.e. taking the trace \begin{align*} \mathrm{div} A = \mathrm{tr}_{12}(\nabla A), \end{align*} where the contraction is to be understood with respect to the metric \(g\). For example, for a 1-form \(\alpha\) and a vector field \(X\) \begin{align*} \mathrm{div}\alpha = g^{ij}(\nabla\alpha)_{ij},\qquad \mathrm{div}X=(\nabla X)^i_{i}. \end{align*} The coordinate expression, therefore, follows directly from the covariant derivative. Again, the covariant divergence of (functions of) the metric vanishes.

[6]:

# vector field
left = Integrate(mf.CovDiv(X), mesh)
right = Integrate(Trace(cov_grad_X), mesh)
print(f"{left:.6f} = {right:.6f}")
assert abs(left - right) < TOL

# 1-form
left = Integrate(mf.CovDiv(alpha), mesh)
right = Integrate(Trace(Inv(metric) * cov_grad_alpha), mesh)
print(f"{left:.6f} = {right:.6f}")
assert abs(left - right) < TOL

# second order tensor field
cov_div_B = Einsum(
    "ij,ijk->k",
    Inv(metric),
    CF((B.Diff(x), B.Diff(y)), dims=(2, 2, 2))
    - Einsum("ijl,lk->ijk", christoffel2, B)
    + Einsum("ilk,jl->ijk", christoffel2, B),
)
left = Integrate(mf.CovDiv(B), mesh).Norm()
right = Integrate(cov_div_B, mesh).Norm()
print(f"{left:.6f} = {right:.6f}")
assert abs(left - right) < TOL

# covariant divergence of (inverse) metric is zero
should_be_zero = Integrate(mf.CovDiv(mf.G), mesh).Norm()
print(f"{should_be_zero:.6f} = 0")
assert should_be_zero < TOL
should_be_zero = Integrate(mf.CovDiv(mf.G_inv), mesh).Norm()
print(f"{should_be_zero:.6f} = 0")
assert should_be_zero < TOL

0.788776 = 0.788776
2.326922 = 2.326922
9.871358 = 9.871358
0.000000 = 0
0.000000 = 0

Analogously to the Euclidean case, there holds the following integration by parts formula for a scalar field \(f \in C^{\infty}(M)\) and a vector field \(X\in\mathfrak{X}(M)\)

\begin{align*} \int_{M} g(\nabla f,X)\,\omega = -\int_{M} f\,\mathrm{div}(X)\,\omega + \int_{\partial M}f\, g(X,n)\,\omega_{\partial M}, \end{align*}

where \(n\) denotes the outer \(g\)-normal vector. That means, \(g(n,n)=1\). Further, for a second order tensor \(A\), we have

\begin{align*} \int_{M} g(\nabla X,A)\,\omega = -\int_{M} g(X,\mathrm{div}(A))\,\omega + \int_{\partial M} g(n\otimes X,A)\,\omega_{\partial M}, \end{align*} where \(g(n\otimes X,A) = g(X, n \llcorner A)\) can be written by contracting \(A\) with \(n\) in its first argument. Note that in the first integration by parts formula \(X\) can be replaced by a 1-form \(\alpha\in \Lambda^1(M)\). Also, \(X\) can be replaced by \(\alpha\) in the second formula, and the second-order tensor \(A\) can have any signature.

More generally, let \(A\in \mathcal{T}^k_l(M)\) and \(B\in\mathcal{T}^r_s(M)\) with \(k+l+1 = r+s\). Then there holds \begin{align*} \int_{M} g(\nabla A,B)\,\omega = -\int_{M} g(A,\mathrm{div}(B))\,\omega + \int_{\partial M} g(n\otimes A,B)\,\omega_{\partial M}. \end{align*}

In coordinates, using that in 2D with \(\hat{t}\) and \(\hat{n}\) the Euclidean tangential and normal vectors \begin{align*} n^i = \frac{\sqrt{\det g}}{\sqrt{g_{\hat{t}\hat{t}}}}g^{ij}\hat{n}_j,\qquad \omega = \sqrt{\det(g)},\qquad \omega_{\partial M} = \sqrt{g_{\hat{t}\hat{t}}} \end{align*} there holds for a scalar \(f\), vector field \(X\) and \((2,0)\)-tensor \(A\)

\begin{align*} \int_M X^i\partial_i f\,\sqrt{\det g}\,dx &= -\int_Mf\,(\partial_iX^i+\Gamma_{ik}^iX^k)\,\sqrt{\det g}\,dx \\ &\quad+ \int_{\partial M} f\,g_{ij}X^i n^j\, \sqrt{g_{tt}}\,ds,\\ \int_M(\partial_iX^j+\Gamma_{ik}^jX^k)g^{ik}A_{kj}\,\sqrt{\det g}\,dx &= -\int_MX^j(\partial_iA_{ij}-\Gamma_{ii}^kA_{kj}-\Gamma_{ij}^kA_{ik})\,\sqrt{\det g}\,dx\\ &\quad +\int_{\partial M}A_{ij}n^i X^j\,\sqrt{g_{tt}}\,ds. \end{align*}

From the above expressions, it is clear that coordinate-based formulas are very error-prone.

[7]:

# scalar with vector field
left = Integrate(mf.InnerProduct(mf.CovDeriv(f), X) * omega_T * dx, mesh)
right = Integrate(
    -mf.InnerProduct(mf.CovDiv(X), f) * omega_T * dx
    + f * mf.InnerProduct(X, mf.normal) * omega_S * ds,
    mesh,
)
print(f"{left:.6f} = {right:.6f}")
assert abs(left - right) < TOL

# vector field with second order tensor field
left = Integrate(mf.InnerProduct(mf.CovDeriv(X), A) * omega_T * dx, mesh)
right = Integrate(
    -mf.InnerProduct(mf.CovDiv(A), X) * omega_T * dx
    + mf.InnerProduct(dg.TensorProduct(mf.normal, X), A) * omega_S * ds,
    mesh,
)
print(f"{left:.6f} = {right:.6f}")
assert abs(left - right) < TOL

0.381579 = 0.381579
1.246865 = 1.246865

Covariant curl and rot in 2D

In two dimensions, the covariant rot and curl operators are defined for scalar fields \(f\), 1-forms \(\alpha\), vector fields \(X\), and \((2,0)\)-tensors \(\sigma\) by \begin{align*} \mathrm{rot}(f)= \varepsilon^{ij}\partial_jf\,\partial_i,\qquad \mathrm{rot}(X) = \varepsilon^{jq}\left(\partial_qX^i+\Gamma^i_{qk}X^k\right)\,\partial_i\otimes \partial_j \end{align*} and \begin{align*} \mathrm{curl}(\alpha)= \varepsilon^{ij}\partial_i\alpha_j,\qquad \mathrm{curl}(\sigma)= \varepsilon^{jk}(\partial_j\sigma_{ik}-\Gamma_{ji}^m\sigma_{mk})\,dx^i. \end{align*} The incompatibility operator of a \((2,0)\)-tensor is defined as two consecutive curl operators \begin{align*} \mathrm{inc}(A) = \mathrm{curl}(\mathrm{curl}(A)),\qquad \mathrm{inc}(A)=\varepsilon^{qi}\varepsilon^{jk}\left(\partial_j\partial_qA_{ik}-\partial_q(\Gamma_{ji}^mA_{mk})-\Gamma_{lq}^l(\partial_jA_{ik}-\Gamma_{ij}^mA_{mk})\right). \end{align*} Above, the covariant Levi-Civita symbols are given by \(\varepsilon^{ij}=1/\sqrt{\det g}\,\hat{\varepsilon}^{ij}\) and \(\varepsilon_{ij}\sqrt{\det g}\,\hat{\varepsilon}_{ij}\), where \(\hat{\varepsilon}_{ij}\) and \(\hat{\varepsilon}_{ij}\) are the standard permuting symbols being 1 or -1 if \((i,j)\) is an even or odd permutation of \((1,2)\) and 0 else. Further, note that the covariant Levi-Civita symbols are not tensors but tensor densities (they get multiplied/divided by the volume form).

We can also apply the rot operator twice in a row \begin{align*} \mathrm{rot}(\mathrm{rot}(f)) = \varepsilon^{jq}\left( \varepsilon^{il}\partial_q\partial_lf-\Gamma_{lq}^q\varepsilon^{il}\partial_lf +\varepsilon^{kl}\Gamma_{qk}^i\partial_lf\right)\,\partial_i\otimes\partial_j. \end{align*}

There holds analogously to the Euclidean case \begin{align*} \mathrm{div}(\mathrm{rot}(f))=0,\qquad \mathrm{curl}(\nabla f) = 0,\qquad \mathrm{curl}(\mathrm{rot}(f)) = -\Delta f, \end{align*} however, these identities are not valid for tensor fields, e.g. \begin{align*} \mathrm{div}(\mathrm{rot}(X))\neq 0,\qquad \mathrm{curl}(\nabla \alpha) \neq 0! \end{align*}

[8]:

# rot
# scalar field
left = Integrate(mf.CovRot(f), mesh).Norm()
right = Integrate(1 / sqrt(Det(metric)) * CF((f.Diff(y), -f.Diff(x))), mesh).Norm()
print(f"{left:.6f} = {right:.6f}")
assert abs(left - right) < TOL

# vector field
cov_rot_X = Einsum(
    "jq,qi->ij",
    1 / sqrt(Det(metric)) * LeviCivitaSymbol(2),
    CF((X.Diff(x), X.Diff(y)), dims=(2, 2)) + Einsum("ikj,k->ij", christoffel2, X),
)
left = Integrate(mf.CovRot(X), mesh).Norm()
right = Integrate(cov_rot_X, mesh).Norm()
print(f"{left:.6f} = {right:.6f}")
assert abs(left - right) < TOL

# 1-form
cov_rot_alpha = Einsum(
    "jq,qi->ij",
    1 / sqrt(Det(metric)) * LeviCivitaSymbol(2),
    dg.GradCF(mf.G_inv[alpha, :], 2)
    + Einsum("ikj,k->ij", christoffel2, mf.G_inv[alpha, :]),
)
left = Integrate(mf.CovRot(alpha), mesh).Norm()
right = Integrate(cov_rot_alpha, mesh).Norm()
print(f"{left:.6f} = {right:.6f}")
assert abs(left - right) < TOL

# curl
# 1-form
cov_curl_alpha = Einsum(
    "ij,ij->",
    1 / sqrt(Det(metric)) * LeviCivitaSymbol(2),
    CF((alpha.Diff(x), alpha.Diff(y)), dims=(2, 2)),
)
left = Integrate(mf.CovCurl(alpha), mesh)
right = Integrate(cov_curl_alpha, mesh)
print(f"{left:.6f} = {right:.6f}")
assert abs(left - right) < TOL

# vector field
cov_curl_X = Einsum(
    "ij,ij->",
    1 / sqrt(Det(metric)) * LeviCivitaSymbol(2),
    dg.GradCF(metric[X, :], 2),
)
left = Integrate(mf.CovCurl(X), mesh)
right = Integrate(cov_curl_X, mesh)
print(f"{left:.6f} = {right:.6f}")
assert abs(left - right) < TOL

# (2,0)-tensor field (general 2-order fields not implemented yet)
cov_curl_C = Einsum(
    "jk,jik->i",
    1 / sqrt(Det(metric)) * LeviCivitaSymbol(2),
    dg.GradCF(C, 2) - Einsum("jim,mk->jik", christoffel2, C),
)
left = Integrate(mf.CovCurl(C), mesh).Norm()
right = Integrate(cov_curl_C, mesh).Norm()
print(f"{left:.6f} = {right:.6f}")
assert abs(left - right) < TOL

1.063312 = 1.063312
10.041697 = 10.041697
3.895205 = 3.895205
-1.145168 = -1.145168
-3.832837 = -3.832837
0.503877 = 0.503877

[9]:

should_be_zero = Integrate(mf.CovDiv(mf.CovRot(f)), mesh)
print(f"div(rot(f)) = {should_be_zero:.6f}")
assert abs(should_be_zero) < TOL

should_be_zero = Integrate(mf.CovCurl(mf.CovDeriv(f)), mesh)
print(f"curl(nabla(f)) = {should_be_zero:.6f}")
assert abs(should_be_zero) < TOL

left = Integrate(mf.CovCurl(mf.CovRot(f)), mesh)
right = -Integrate(mf.Trace(mf.CovHesse(f)), mesh)
print(f"curl(rot(f)) = -Delta f: {left:.6f} = {right:.6f}")
assert abs(left - right) < TOL

# invoking tensor fields
print(
    f"curl(nabla(alpha)) = {Integrate(mf.CovCurl(mf.CovDeriv(alpha)), mesh).Norm():.6f}"
)
print(f"div(rot(X)) = {Integrate(mf.CovDiv(mf.CovRot(X)), mesh).Norm():.6f}")

div(rot(f)) = 0.000000
curl(nabla(f)) = -0.000000
curl(rot(f)) = -Delta f: -1.833333 = -1.833333
curl(nabla(alpha)) = 2.824292
div(rot(X)) = 18.249655

With \(t\) denoting the \(g\)-tangential vector, in coordinates \(t= \frac{1}{\sqrt{g_{\hat{t}\hat{t}}}}\hat{t}^i\,\partial_i\), there holds the following integration by parts formulae \begin{align*} \int_M\mathrm{curl}(X)\,f\,\omega &= \int_M g(\mathrm{rot}(f),X)\,\omega + \int_{\partial M} g(X,t)\,f\,\omega_{\partial M},\\ \int_Mg(\mathrm{curl}(A),X)\,\omega &= \int_M g(A,\mathrm{rot}(X))\,\omega + \int_{\partial M}g(A,X\otimes t)\,\omega_{\partial M},\\ \int_M \mathrm{inc}(A)\,f\,\omega&=\int_{M}g(\mathrm{curl}(A),\mathrm{rot}(f))\,\omega +\int_{\partial M}g(\mathrm{curl}(A),t)\,f\,\omega_{\partial M}\\ &= \int_Mg(A,\mathrm{rot}(\mathrm{rot}(f)))\,\omega + \int_{\partial M}\left( g(A,\mathrm{rot}(f)\otimes t)+ g(\mathrm{curl}(A),t)\,f\right)\omega_{\partial M}. \end{align*}

[10]:

# curl-rot integration-by-parts between vector and scalar field
left = Integrate(mf.InnerProduct(mf.CovCurl(X), f) * omega_T * dx, mesh)
right = Integrate(
    mf.InnerProduct(X, mf.CovRot(f)) * omega_T * dx
    + f * mf.InnerProduct(X, mf.tangent) * omega_S * ds,
    mesh,
)
print(f"{left:.6f} = {right:.6f}")
assert abs(left - right) < TOL

# curl-rot integration-by-parts between (2,0)-tensor and vector field
left = Integrate(mf.InnerProduct(mf.CovCurl(C), X) * omega_T * dx, mesh)
right = Integrate(
    mf.InnerProduct(C, mf.CovRot(X)) * omega_T * dx
    + mf.InnerProduct(C, dg.TensorProduct(X, mf.tangent)) * omega_S * ds,
    mesh,
)
print(f"{left:.6f} = {right:.6f}")
assert abs(left - right) < TOL

# inc - rot/rot integration-by-parts between (2,0)-tensor and scalar field
left = Integrate(mf.InnerProduct(mf.CovInc(C, True), f) * omega_T * dx, mesh)
mid = Integrate(
    mf.InnerProduct(mf.CovCurl(C), mf.CovRot(f)) * omega_T * dx
    + f
    * mf.InnerProduct(mf.CovCurl(C), mf.tangent)
    * omega_S
    * dx(element_boundary=True),
    mesh,
)
right = Integrate(
    mf.InnerProduct(C, mf.CovRot(mf.CovRot(f))) * omega_T * dx
    + (
        mf.InnerProduct(C, dg.TensorProduct(mf.CovRot(f), mf.tangent))
        + f * mf.InnerProduct(mf.CovCurl(C), mf.tangent)
    )
    * omega_S
    * dx(element_boundary=True),
    mesh,
)
print(f"{left:.6f} = {mid:.6f} = {right:.6f}")
assert abs(left - mid) < TOL
assert abs(left - right) < TOL

-0.562212 = -0.562212
0.610140 = 0.610140
-1.290332 = -1.290332 = -1.290332

With our definition that the normal vector \(n\) points outward elements, there holds the following relation with the tangent vector \(t\) \begin{align*} g(\mathrm{rot}(f),n) = g(\nabla f,t). \end{align*}

[11]:

left = Integrate(
    mf.InnerProduct(mf.CovRot(f), mf.normal) * dx(element_boundary=True), mesh
)
right = Integrate(
    mf.InnerProduct(mf.CovDeriv(f), mf.tangent) * dx(element_boundary=True), mesh
)
print(f"{left:.6f} = {right:.6f}")
assert abs(left - right) < TOL

-0.053979 = -0.053979

Covariant curl in 3D

Next, we consider 3-dimensional Riemannian manifolds. First, define a metric and some fields.

[12]:

metric = dg.WarpedProduct().metric
mesh = Mesh(unit_cube.GenerateMesh(maxh=0.1))
mf = dg.RiemannianManifold(metric)
omega_T = mf.VolumeForm(VOL)
omega_S = mf.VolumeForm(BND)

f = dg.ScalarField(CF(x**2 * y - 0.1 * y * x + 0.2 * y * z**2), dim=3)
X = dg.VectorField(CF((10 * x * y**3 - x**2 * z, y**4 * x - y * z, z * x * y)))
Y = dg.VectorField(
    CF(
        (
            y**2 * x * z - 0.1 * y * x**2,
            10 * x * y**3 + x**2 - y * z,
            sin(z) + cos(x * y),
        )
    )
)

alpha = dg.OneForm(CF((sin(x * y), cos(x) * y**2, 2 * x * y * z)))
beta = dg.OneForm(
    CF(
        (
            0.2 * y * x**2 + 0.37 * x**3 * z,
            2 * x**2 * y**2 + 0.1 * x**2 - 1.34 * y * z,
            cos(x * y) + 0.1 * z,
        )
    )
)

B = dg.TensorProduct(alpha, X)

The Covariant curl is defined with the covariant Levi-Civita symbols \(\varepsilon^{ijk}=1/\sqrt{\det g}\,\hat{\varepsilon}^{ijk}\) and \(\varepsilon_{ijk}=\sqrt{\det g}\,\hat{\varepsilon}_{ijk}\) for vector fields \(X\in\mathfrak{X}(M)\) and 1-forms \(\alpha\in\Lambda^1(M)\) \begin{align*} \mathrm{curl}(X) = \frac{1}{\sqrt{\det g}}\varepsilon^{ijk}\partial{j}X_k\,\partial_i,\qquad \mathrm{curl}(\alpha) = \frac{1}{\sqrt{\det g}}\varepsilon^{ijk}\partial_{j}\alpha_k\,\partial_i, \end{align*} where we lowered the index for \(X\).

There holds the integration by parts formula \begin{align*} \int_M g(\mathrm{curl}(X),\alpha)\,\omega = \int_M g(X,\mathrm{curl}(\alpha))\,\omega - \int_{\partial M}g(X\times n,\alpha)\,\omega_{\partial M}, \end{align*} where \(n\) is again the \(g\)-normal vector and the cross product is defined by \(X\times Y = \omega(X,Y,\cdot)^\sharp\), or in coordinates \begin{align*} X\times Y = \frac{1}{\sqrt{\det g}}\varepsilon^{ijk}X_jY_k\,\partial_i. \end{align*} We can also use 1-forms and mix vector fields with 1-forms in the above formula.

[13]:

with TaskManager():
    # two vector fields
    left = Integrate(
        mf.InnerProduct(mf.CovCurl(X), Y) * omega_T * dx(bonus_intorder=1), mesh
    )
    right = Integrate(
        mf.InnerProduct(X, mf.CovCurl(Y)) * omega_T * dx(bonus_intorder=1)
        - mf.InnerProduct(mf.Cross(X, mf.normal), Y) * omega_S * ds(bonus_intorder=1),
        mesh,
    )
    print(f"{left:.6f} = {right:.6f}")
    assert abs(left - right) < TOL

    # vector field with 1-form
    left = Integrate(
        mf.InnerProduct(mf.CovCurl(X), alpha) * omega_T * dx(bonus_intorder=1), mesh
    )
    right = Integrate(
        mf.InnerProduct(X, mf.CovCurl(alpha)) * omega_T * dx(bonus_intorder=1)
        - mf.InnerProduct(mf.Cross(X, mf.normal), alpha)
        * omega_S
        * ds(bonus_intorder=1),
        mesh,
    )
    print(f"{left:.6f} = {right:.6f}")
    assert abs(left - right) < TOL

    # 1-form fields
    left = Integrate(
        mf.InnerProduct(mf.CovCurl(alpha), beta) * omega_T * dx(bonus_intorder=1), mesh
    )
    right = Integrate(
        mf.InnerProduct(alpha, mf.CovCurl(beta)) * omega_T * dx(bonus_intorder=1)
        - mf.InnerProduct(mf.Cross(alpha, mf.normal), beta)
        * omega_S
        * ds(bonus_intorder=1),
        mesh,
    )
    print(f"{left:.6f} = {right:.6f}")
    assert abs(left - right) < TOL

78.498239 = 78.498239
-5.844097 = -5.844097
-0.400652 = -0.400652

Note, that we used the bonus-intorder flag to increase the order of numerical integration to integrate the terms sufficiently accurate.

Hessian and Laplace-Beltrami

For a scalar field \(f\in C^{\infty}(M)\) the Riemannian Hessian is defined by \(\nabla^2 f = \nabla df\in \mathcal{T}^2_0(M)\). In coordinates

\begin{align*} \mathrm{Hess}f = \left(\partial^2_{i,j}f-\Gamma^k_{ij}\partial_k f\right)\,dx^i\otimes dx^j. \end{align*} The Hessian is symmetric by the symmetry of Christoffel symbols \(\Gamma_{ij}^k=\Gamma_{ji}^k\).

The Laplace-Beltrami operator of a scalar field \(f\in C^{\infty}(M)\) is defined in the usual way \(\Delta f = \mathrm{div}(\nabla f)\), in coordinates \begin{align*} \Delta f = \mathrm{tr}(\nabla^2 f) = \partial^2_{i,i}f-\Gamma^k_{ii}\partial_k f = \frac{1}{\sqrt{\det g}}\mathrm{div}(\sqrt{\det g} (\nabla f)^\sharp). \end{align*}

There hold the following integration by parts formulae \begin{align*} &\int_M g(\mathrm{div}A,\nabla f)\,\omega = - \int_Mg(A,\mathrm{Hess}f)\,\omega + \int_{\partial M} g(A, n\otimes \nabla f)\omega_{\partial M},\\ &\int_M \Delta f\,g\,\omega = -\int_{M} g(\nabla f,\nabla g)\,\omega + \int_{\partial M}(\nabla f)(n)\,g\,\omega_{\partial M}. \end{align*} Note that \((\nabla f)(n) = g(\nabla f, n)\).

[14]:

with TaskManager():
    left = Integrate(mf.InnerProduct(mf.CovDiv(B), mf.CovDeriv(f)) * omega_T * dx, mesh)
    right = Integrate(
        -mf.InnerProduct(B, mf.CovHesse(f)) * omega_T * dx
        + mf.InnerProduct(B, dg.TensorProduct(mf.normal, mf.CovDeriv(f)))
        * omega_S
        * dx(element_boundary=True),
        mesh,
    )
    print(f"{left:.6f} = {right:.6f}")
    assert abs(left - right) < TOL

16.867421 = 16.867421

Note that we used dx(element-boundary=True) instead of ds above for the right-hand side. Otherwise, the full covariant derivative would not be available.

[15]:

with TaskManager():
    left = Integrate(mf.InnerProduct(mf.Trace(mf.CovHesse(f)), g) * omega_T * dx, mesh)
    right = Integrate(
        -mf.InnerProduct(mf.CovDeriv(f), mf.CovDeriv(g)) * omega_T * dx
        + mf.InnerProduct(mf.CovDeriv(f), mf.normal)
        * g
        * omega_S
        * dx(element_boundary=True),
        mesh,
    )
    print(f"{left:.6f} = {right:.6f}")
    assert abs(left - right) < TOL

-0.202355 = -0.202355

Fancy identities

2D

For a \((2,0)\)-tensor \(\sigma\in \mathcal{T}^{2}_0(M)\) there holds \begin{align*} \mathrm{inc}(\sigma) = -\mathrm{div}\mathrm{div}(\mathbb{S}\sigma) = -\mathrm{div}\mathrm{div}(\sigma) + \Delta \mathrm{Tr}(\sigma), \end{align*} where \(\mathbb{S}\sigma = \sigma-\mathrm{Tr}(\sigma)\,g\).

[16]:

metric = dg.CigarSoliton().metric
mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))
mf = dg.RiemannianManifold(metric)

alpha = dg.OneForm(CF((sin(x * y), cos(x) * y**2)))
beta = dg.OneForm(
    CF((0.2 * y * x**2 + 0.37 * x**3, 2 * x**2 * y**2 + 0.1 * x**2 - 1.34 * y))
)

C = dg.TensorProduct(alpha, beta)


lhs = Integrate(mf.CovInc(C, True), mesh)
mid = Integrate(
    -mf.CovDiv(mf.CovDiv(dg.TensorField(C - mf.Trace(C) * mf.G, "11"))), mesh
)
rhs = Integrate(-mf.CovDiv(mf.CovDiv(C)) + mf.Trace(mf.CovHesse((mf.Trace(C)))), mesh)

print(f"{lhs:.6f} = {mid:.6f} = {rhs:.6f}")
assert abs(lhs - mid) < TOL
assert abs(mid - rhs) < TOL

-4.573159 = -4.573159 = -4.573159

[ ]: