Alternating k-forms, wedge product, Hodge star, and exterior (co-) derivative

In this notebook, we introduce important objects and identities of exterior calculus.

[1]:

from ngsolve import *
from netgen.occ import unit_square, unit_cube
import ngsdiffgeo as dg

TOL = 1e-8

def l2_error(a, b, mesh):
    return sqrt(Integrate(InnerProduct(a - b, a - b), mesh))

def l2_norm(a, mesh):
    return sqrt(Integrate(InnerProduct(a, a), mesh))

3D Euclidean case

We focus on the more interesting case of \(N=3\). First we consider the Euclidean manifold, \(\Omega\subset \mathbb{R}^3\), with the usual Euclidean inner products. Although the wedge product and exterior derivative are independent of the metric, for the Hodge star and exterior co-derivative this is essential. Later, we consider a non-Euclidean Riemannian manifold with a metric tensor \(g\neq I\).

We define the 1-forms \(dx,dy,dz\in \Lambda^1(\Omega)\), 2-forms \(dx\wedge dy, dx\wedge dz, dy\wedge dz\in \Lambda^2(\Omega)\), and 3-form \(dx\wedge dy\wedge dz\in \Lambda^3(\Omega)\). The wedge product anti-symmetrizes the dyadic product of the forms. For example there holds

\[\begin{align*} \alpha\wedge \beta = \alpha \otimes \beta - \beta\otimes \alpha,\qquad \alpha,\beta\in\Lambda^1(\Omega). \end{align*}\]

For the wedge product, there holds the associative rule \((\alpha\wedge\beta)\wedge \gamma = \alpha\wedge(\beta\wedge\gamma)\).

[2]:

mesh3 = Mesh(unit_cube.GenerateMesh(maxh=2))

dx = dg.OneForm(CF((1, 0, 0)))
dy = dg.OneForm(CF((0, 1, 0)))
dz = dg.OneForm(CF((0, 0, 1)))

dx_dy = dg.Wedge(dx, dy)
dx_dz = dg.Wedge(dx, dz)
dy_dz = dg.Wedge(dy, dz)

dx_dy2 = dg.TwoForm(CF((0, 1, 0, -1, 0, 0, 0, 0, 0), dims=(3, 3)), dim=3)
dx_dz2 = dg.TwoForm(CF((0, 0, 1, 0, 0, 0, -1, 0, 0), dims=(3, 3)), dim=3)
dy_dz2 = dg.TwoForm(CF((0, 0, 0, 0, 0, 1, 0, -1, 0), dims=(3, 3)), dim=3)

print(f"L2 error dx_dy - dx_dy2 = {l2_error(dx_dy, dx_dy2, mesh3):.6f}")
assert l2_error(dx_dy, dx_dy2, mesh3) < TOL
print(f"L2 error dx_dz - dx_dz2 = {l2_error(dx_dz, dx_dz2, mesh3):.6f}")
assert l2_error(dx_dz, dx_dz2, mesh3) < TOL
print(f"L2 error dy_dz - dy_dz2 = {l2_error(dy_dz, dy_dz2, mesh3):.6f}")
assert l2_error(dy_dz, dy_dz2, mesh3) < TOL

dx_dy_dz = dg.Wedge(dx_dy, dz)
dx_dy_dz_as = dg.Wedge(dx, dy_dz)
dx_dy_dz2 = dg.ThreeForm(
    CF(
        (
            0,
            0,
            0,
            0,
            0,
            1,
            0,
            -1,
            0,
            0,
            0,
            -1,
            0,
            0,
            0,
            1,
            0,
            0,
            0,
            1,
            0,
            -1,
            0,
            0,
            0,
            0,
            0,
        ),
        dims=(3, 3, 3),
    ),
    dim=3,
)
print(f"L2 error dx_dy_dz - dx_dy_dz2 = {l2_error(dx_dy_dz, dx_dy_dz2, mesh3):.6f}")
assert l2_error(dx_dy_dz, dx_dy_dz2, mesh3) < TOL
print(f"L2 error dx_dy_dz - dx_dy_dz_as = {l2_error(dx_dy_dz, dx_dy_dz_as, mesh3):.6f}")
assert l2_error(dx_dy_dz, dx_dy_dz_as, mesh3) < TOL

L2 error dx_dy - dx_dy2 = 0.000000
L2 error dx_dz - dx_dz2 = 0.000000
L2 error dy_dz - dy_dz2 = 0.000000
L2 error dx_dy_dz - dx_dy_dz2 = 0.000000
L2 error dx_dy_dz - dx_dy_dz_as = 0.000000

Further, the wedge product is anti-commutative depending on the degree of the involved forms

\[\begin{align*} \alpha\wedge\beta = (-1)^{lk}\beta\wedge\alpha,\qquad \alpha\in\Lambda^l(\Omega), \beta\in\Lambda^k(\Omega), \end{align*}\]

such that \(\alpha\wedge\alpha=0\) for all \(\alpha\in\Lambda^k(\Omega)\).

[3]:

dy_dx = dg.Wedge(dy, dx)
print(f"L2 error dy^dx + dx^dy = {l2_norm(dy_dx + dx_dy, mesh3):.6f}")
assert l2_norm(dy_dx + dx_dy, mesh3) < TOL

dz_dx_dy = dg.Wedge(dz, dx_dy)
print(f"L2 error dz^(dx^dy) - dx^dy^dz = {l2_norm(dz_dx_dy - dx_dy_dz, mesh3):.6f}")
assert l2_norm(dz_dx_dy - dx_dy_dz, mesh3) < TOL

alpha = dg.OneForm(CF((0.3 * x * y, z**2, -0.1 * x)))
alpha_wedge_alpha = dg.Wedge(alpha, alpha)
print(f"L2 norm alpha^alpha = {l2_norm(alpha_wedge_alpha, mesh3):.6f}")
assert l2_norm(alpha_wedge_alpha, mesh3) < TOL

L2 error dy^dx + dx^dy = 0.000000
L2 error dz^(dx^dy) - dx^dy^dz = 0.000000
L2 norm alpha^alpha = 0.000000

Let \(X\in\mathfrak{X}(\Omega)\) be a vector field. The contraction (also denoted as interior product) \(\iota_X:\Lambda^k(\Omega)\to\Lambda^{k-1}(\Omega)\), sometimes denoted by \(X\lrcorner\), is defined by

\[\begin{align*} (\iota_X\alpha)(Y_1,\dots,Y_{k-1}) = \alpha(X,Y_1,\dots, Y_{k-1}), \qquad\alpha\in\Lambda^k(\Omega),\, Y_i\in \mathfrak{X}(\Omega). \end{align*}\]

There holds the following identity for the wedge product

\[\begin{align*} \iota_X(\beta\wedge\gamma) = (\iota_X\beta)\wedge \gamma + (-1)^k\beta\wedge(\iota_X\gamma),\qquad\beta\in\Lambda^k(\Omega),\gamma\in\Lambda^l(\Omega). \end{align*}\]

[4]:

mf = dg.RiemannianManifold(Id(3))

X = dg.VectorField(CF((0.3 * y * z - x * y, x**2 * z + 0.34 * y**3, -x * y * z)))
alpha = dg.OneForm(CF((0.3 * x * y, z**2, -0.1 * x)))
beta = dg.TwoForm(
    CF(
        (
            0,
            y * z - 0.1 * x**2,
            0,
            -(y * z - 0.1 * x**2),
            0,
            x * z,
            0,
            -x * z,
            0,
        ),
        dims=(3, 3),
    )
)

iota_X_beta_gamma = mf.Contraction(X, dg.Wedge(alpha, beta))
iota_beta_iota_gamma = dg.Wedge(mf.Contraction(X, alpha), beta) - dg.Wedge(
    alpha, mf.Contraction(X, beta)
)
print(
    f"L2 error i_X(alpha^beta) - i_X(alpha)^beta + alpha^i_X(beta) = {l2_error(iota_X_beta_gamma, iota_beta_iota_gamma, mesh3):.6f}"
)
assert l2_error(iota_X_beta_gamma, iota_beta_iota_gamma, mesh3) < TOL

L2 error i_X(alpha^beta) - i_X(alpha)^beta + alpha^i_X(beta) = 0.000000

The Hodge star operator \(\star:\Lambda^{k}(\Omega)\to \Lambda^{N-k}(\Omega)\) is defined by the identity

\[\begin{align*} \alpha\wedge (\star\beta) = \frac{1}{k!}\langle \alpha, \beta\rangle\,\omega,\qquad \alpha,\beta\in \Lambda^{k}(\Omega), \end{align*}\]

where \(\omega\) is the volume form (\(dx\wedge dy\wedge dz\) in the Euclidean case) and \(\langle \alpha, \beta\rangle\) the inner product of \(k\)-tensors. Note that the normalization factor \(1/k!\) is necessary as in the inner product we sum over all indices. In the literature the inner product of \(k\)-forms, sometimes denoted by \(g^{1}(\alpha,\beta)\) is defined by

\[\begin{align*} g^{-1}(\alpha, \beta)=\frac{1}{k!} g^{i_1 j_1}\cdots g^{i_k j_k}\, \alpha_{i_1\dots i_k}\, \beta_{j_1\dots j_k},\qquad \alpha,\beta\in \Lambda^{k}(\Omega), \end{align*}\]

without the normalization factor, so there holds \(\langle\alpha,\beta\rangle = k!\,g^{-1}(\alpha,\beta)\) for \(\alpha,\beta\in\Lambda^k(\Omega)\). Applying the Hodge star twice yields the identity up to a sign

\[\begin{align*} \star\star \alpha = (-1)^{k(N-k)}\alpha, \alpha\in\Lambda^k(\Omega). \end{align*}\]

[5]:

# Euclidean manifold, metric tensor g=Id
mf = dg.RiemannianManifold(Id(3))

one = dg.ScalarField(CF(1), dim=3)

star_one = mf.star(one)
star_dx = mf.star(dx)
star_dy = mf.star(dy)
star_dz = mf.star(dz)
star_dx_dy = mf.star(dx_dy)
star_dx_dz = mf.star(dx_dz)
star_dy_dz = mf.star(dy_dz)
star_dx_dy_dz = mf.star(dx_dy_dz)

print(f"L2 error star(1) - dx^dy^dz = {l2_error(star_one, dx_dy_dz, mesh3):.6f}")
assert l2_error(star_one, dx_dy_dz, mesh3) < TOL
print(f"L2 error star(dx) - dy^dz = {l2_error(star_dx, dy_dz, mesh3):.6f}")
assert l2_error(star_dx, dy_dz, mesh3) < TOL
print(f"L2 error star(dy) - -dz^dx = {l2_error(star_dy, -dx_dz, mesh3):.6f}")
assert l2_error(star_dy, -dx_dz, mesh3) < TOL
print(f"L2 error star(dz) - dx^dy = {l2_error(star_dz, dx_dy, mesh3):.6f}")
assert l2_error(star_dz, dx_dy, mesh3) < TOL
print(f"L2 error star(dx^dy) - dz = {l2_error(star_dx_dy, dz, mesh3):.6f}")
assert l2_error(star_dx_dy, dz, mesh3) < TOL
print(f"L2 error star(dx^dz) - -dy = {l2_error(star_dx_dz, -dy, mesh3):.6f}")
assert l2_error(star_dx_dz, -dy, mesh3) < TOL
print(f"L2 error star(dy^dz) - dx = {l2_error(star_dy_dz, dx, mesh3):.6f}")
assert l2_error(star_dy_dz, dx, mesh3) < TOL
print(f"L2 error star(dx^dy^dz) - 1 = {l2_error(star_dx_dy_dz, one, mesh3):.6f}")
assert l2_error(star_dx_dy_dz, one, mesh3) < TOL

alpha = dg.OneForm(CF((0.3 * x * y, z**2, -0.1 * x)))
beta = dg.TwoForm(
    CF(
        (
            0,
            y * z - 0.1 * x**2,
            0,
            -(y * z - 0.1 * x**2),
            0,
            x * z,
            0,
            -x * z,
            0,
        ),
        dims=(3, 3),
    )
)
gamma = (0.1 + x * y - 0.2 * x * y * z) * dx_dy_dz

star_star_alpha = mf.star(mf.star(alpha))
expected_star_star_alpha = alpha
star_star_beta = mf.star(mf.star(beta))
expected_star_star_beta = beta
star_star_gamma = mf.star(mf.star(gamma))
expected_star_star_gamma = gamma

print(
    f"L2 error star(star(alpha)) - alpha = {l2_error(star_star_alpha, expected_star_star_alpha, mesh3):.6f}"
)
assert l2_error(star_star_alpha, expected_star_star_alpha, mesh3) < TOL
print(
    f"L2 error star(star(beta)) - beta = {l2_error(star_star_beta, expected_star_star_beta, mesh3):.6f}",
)
assert l2_error(star_star_beta, expected_star_star_beta, mesh3) < TOL
print(
    f"L2 error star(star(gamma)) - gamma = {l2_error(star_star_gamma, expected_star_star_gamma, mesh3):.6f}",
)
assert l2_error(star_star_gamma, expected_star_star_gamma, mesh3) < TOL

L2 error star(1) - dx^dy^dz = 0.000000
L2 error star(dx) - dy^dz = 0.000000
L2 error star(dy) - -dz^dx = 0.000000
L2 error star(dz) - dx^dy = 0.000000
L2 error star(dx^dy) - dz = 0.000000
L2 error star(dx^dz) - -dy = 0.000000
L2 error star(dy^dz) - dx = 0.000000
L2 error star(dx^dy^dz) - 1 = 0.000000
L2 error star(star(alpha)) - alpha = 0.000000
L2 error star(star(beta)) - beta = 0.000000
L2 error star(star(gamma)) - gamma = 0.000000

The exterior derivative \(d:\Lambda^k(\Omega)\to\Lambda^{k+1}(\Omega)\) is defined by its action on vector fields \(X_0,\dots,X_k\in\mathfrak{X}(\Omega)\)

\[\begin{align*} (d\alpha)(X_0,\dots,X_k) = \sum_{i=0}^k (-1)^i X_i\big(\alpha(X_0,\dots,\hat{X_i},\dots,X_k)\big) + \sum_{0\leq i<j\leq k} (-1)^{i+j} \alpha([X_i,X_j],X_0,\dots,\hat{X_i},\dots,\hat{X_j},\dots,X_k), \end{align*}\]

which in coordinates reduces to the alternating sum of partial derivatives of the coefficients. Here, \([X,Y]\) denotes the Lie-bracket (or commutator), which is defined by the identity \([X,Y]f = X(Y(f))-Y(X(f))\) for all functions \(f\). Further, \(\hat X\) means that this argument gets skipped.

Applying the exterior derivative twice yields zero, \(dd \alpha = 0\) for \(\alpha\in\Lambda^k(\Omega)\), i.e., \(d\) is nilpotent. Further, there holds the following Leibnitz rule with the wedge product

\[\begin{align*} d(\alpha\wedge\beta) = d\alpha\wedge\beta + (-1)^k\alpha\wedge d\beta,\quad\alpha\in\Lambda^k(\Omega),\beta\in\Lambda^l(\Omega). \end{align*}\]

For 1-forms there also holds \(d \alpha = \nabla \alpha - (\nabla\alpha)^T\) with \(\nabla\) denoting the Levi-Civita connection (covariant derivative) as the involved Christoffel symbols cancel out due to their symmetry.

[6]:

g = dg.ScalarField(x * y * z, dim=3)
alpha = dg.OneForm(CF((0.3 * x * y, z**2, -0.1 * x)))
beta = dg.OneForm(CF((0.3 * z * y, x * z**2, y**2)))

dg_g = dg.d(g)
expected_dg_g = dg.OneForm(CF((y * z, x * z, x * y)))
print(f"|d g - expected||_L2: {l2_norm(dg_g - expected_dg_g, mesh3):.6f}")
assert l2_norm(dg_g - expected_dg_g, mesh3) < TOL
ddg_g = dg.d(dg_g)
print(f"d^2 g = {l2_norm(ddg_g, mesh3):.6f}")
assert l2_norm(ddg_g, mesh3) < TOL

dalpha = dg.d(alpha)
expected_dalpha = mf.CovDeriv(alpha) - mf.CovDeriv(alpha).trans
print(f"||d alpha - expected||_L2: {l2_norm(dalpha - expected_dalpha, mesh3):.6f}")
assert l2_norm(dalpha - expected_dalpha, mesh3) < TOL
ddalpha = dg.d(dalpha)
print(f"d^2 alpha = {l2_norm(ddalpha, mesh3):.6f}")
assert l2_norm(ddalpha, mesh3) < TOL

d_a_w_b = dg.d(dg.Wedge(alpha, beta))
expected_d_a_w_b = dg.Wedge(dalpha, beta) - dg.Wedge(alpha, dg.d(beta))
print(
    f"||d(alpha^beta) - expected||_L2: {l2_norm(d_a_w_b - expected_d_a_w_b, mesh3):.6f}"
)
assert l2_norm(d_a_w_b - expected_d_a_w_b, mesh3) < TOL

|d g - expected||_L2: 0.000000
d^2 g = 0.000000
||d alpha - expected||_L2: 0.000000
d^2 alpha = 0.000000
||d(alpha^beta) - expected||_L2: 0.000000

For a vector field \(X\in\mathfrak{X}(\Omega)\) the covariant divergence and curl (in 3D) can be written in terms of the exterior derivative, Hodge star, and the musical isomorphisms relating 1-forms with vector fields

\[\begin{align*} \mathrm{div}\, X = \star d\star (X^\flat),\qquad \mathrm{curl}\,X = (\star d(X^\flat))^\sharp. \end{align*}\]

[7]:

mf = dg.RiemannianManifold(Id(3))
X = dg.VectorField(CF((0.3 * y * z - x * y, x**2 * z + 0.34 * y**3, -x * y * z)))

div_X = mf.star(dg.d(mf.star(mf.Lower(X))))
div_X_expected = mf.CovDiv(X)
print(f"||div X - expected||_L2: {l2_norm(div_X - div_X_expected, mesh3):.6f}")
assert l2_norm(div_X - div_X_expected, mesh3) < TOL

curl_X = mf.Raise(mf.star(dg.d(mf.Lower(X))))
curl_X_expected = mf.CovCurl(X)
print(f"||curl X - expected||_L2: {l2_norm(curl_X - curl_X_expected, mesh3):.6f}")
assert l2_norm(curl_X - curl_X_expected, mesh3) < TOL

||div X - expected||_L2: 0.000000
||curl X - expected||_L2: 0.000000

Further, there holds the Stokes theorem for \(\alpha\in \Lambda^{N-1}(\Omega)\)

\[\begin{align*} \int_{\Omega}d\alpha = \int_{\partial\Omega}\iota^*\alpha, \end{align*}\]

where \(\iota^*\) is the pullback of the inclusion \(\iota:\partial\Omega\to\Omega\). We use that \(d\alpha = \star(d\alpha)\omega\) and \(\iota^*\alpha = \star_\partial(\iota^*\alpha)\omega_{\partial}\) for \(\alpha\in \Lambda^{N-1}(\Omega)\).

[8]:

import ngsolve

mf = dg.RiemannianManifold(Id(3))
beta = dg.TwoForm(
    CF(
        (
            0,
            y * z - 0.1 * x**2,
            0,
            -(y * z - 0.1 * x**2),
            0,
            0.2 * x * z,
            0,
            -0.2 * x * z,
            0,
        ),
        dims=(3, 3),
    )
)


left = Integrate(mf.star(dg.d(beta)) * ngsolve.dx(bonus_intorder=5), mesh3)
right = Integrate(
    mf.star(beta, BND) * ngsolve.dx(bonus_intorder=5, element_boundary=True), mesh3
)

print(f"Left integral: {left:.6f}, right integral: {right:.6f}")
assert abs(left - right) < TOL

Left integral: 0.600000, right integral: 0.600000

The exterior coderivative \(\delta:\Lambda^k(\Omega)\to\Lambda^{k-1}(\Omega)\) is defined as the adjoint of the exterior derivative with respect to the \(L^2\) inner product, i.e., for all \(\alpha\in\Lambda^k(\Omega)\) and \(\beta\in \Lambda^{k+1}(\Omega)\) with compact support there holds

\[\begin{align*} \int_{\Omega}\frac{1}{(k+1)!}\langle d\alpha,\beta\rangle\,\omega = \int_{\Omega}\frac{1}{k!}\langle \alpha,\delta\beta\rangle\,\omega, \end{align*}\]

where \(\omega\) and \(\langle\cdot,\cdot\rangle\) are again the volume form and the inner product for \(k\)-tensor fields. There holds the identity \(\delta = (-1)^{N(k+1)+1}\,\star d \star\). From it, there follows that \(\delta\) is also nilpotent, \(\delta\delta=0\).

Using Stokes’ theorem together with yields the following integration-by-parts formula that holds for all \(\alpha\in\Lambda^k(\Omega)\) and \(\beta\in \Lambda^{k+1}(\Omega)\)

\[\begin{align*} \int_{\Omega}\frac{1}{(k+1)!}\langle d\alpha,\beta\rangle\,\omega = \int_{\Omega}\frac{1}{k!}\langle \alpha,\delta\beta\rangle\,\omega + \int_{\partial \Omega}\frac{1}{k!}\langle \alpha,\iota_n\beta\rangle\,\omega_{\partial}, \end{align*}\]

where \(\omega_{\partial}\) is the volume form on the boundary, \(n\) the outer unit normal vector, and for a vector field \(X\in\mathfrak{X}(\Omega)\) the contraction is \(\iota_X\). Further, the inner product in the boundary term has to be understood with respect to the tangential metric \(g_F\).

[9]:

import ngsolve

mf = dg.RiemannianManifold(Id(3))

f = dg.ScalarField(200 * x * (1 - x) * y * (1 - y) * z * (1 - z), dim=3)
alpha = dg.OneForm(CF((0.3 * x * y, z**2, -0.1 * x)))

left = Integrate(
    mf.InnerProduct(alpha, dg.d(f)) * mf.VolumeForm(VOL) * ngsolve.dx(bonus_intorder=5),
    mesh3,
)
right = Integrate(
    mf.InnerProduct(f, mf.delta(alpha))
    * mf.VolumeForm(VOL)
    * ngsolve.dx(bonus_intorder=5),
    mesh3,
)

print(f"Compact support: left: {left:.6f}, right: {right:.6f}")
assert abs(left - right) < TOL

f = dg.ScalarField(20 * x * y * (1 - y) * (1 - z), dim=3)
alpha = dg.OneForm(CF((0.3 * x * y, z**2, -0.1 * x)))

left = Integrate(
    mf.InnerProduct(alpha, dg.d(f)) * mf.VolumeForm(VOL) * ngsolve.dx(bonus_intorder=5),
    mesh3,
)
right = Integrate(
    mf.InnerProduct(f, mf.delta(alpha))
    * mf.VolumeForm(VOL)
    * ngsolve.dx(bonus_intorder=5),
    mesh3,
) + Integrate(
    f
    * mf.Contraction(mf.normal, alpha)
    * mf.VolumeForm(BND)
    * ngsolve.dx(element_boundary=True, bonus_intorder=5),
    mesh3,
)
print(f"Integration by parts: left: {left:.6f}, right: {right:.6f}")
assert abs(left - right) < TOL

alpha = dg.OneForm(CF((0.3 * x * y, z**2, -0.1 * x)))
beta = dg.TwoForm(
    CF(
        (
            0,
            y * z - 0.1 * x**2,
            0,
            -(y * z - 0.1 * x**2),
            0,
            x * z,
            0,
            -x * z,
            0,
        ),
        dims=(3, 3),
    )
)
left = Integrate(
    0.5
    * mf.InnerProduct(beta, dg.d(alpha))
    * mf.VolumeForm(VOL)
    * ngsolve.dx(bonus_intorder=5),
    mesh3,
)
right = Integrate(
    mf.InnerProduct(alpha, mf.delta(beta))
    * mf.VolumeForm(VOL)
    * ngsolve.dx(bonus_intorder=5),
    mesh3,
) + Integrate(
    mf.InnerProduct(alpha, mf.Contraction(mf.normal, beta), BND)
    * mf.VolumeForm(BND)
    * ngsolve.dx(element_boundary=True, bonus_intorder=5),
    mesh3,
)
print(f"Integration by parts: left: {left:.6f}, right: {right:.6f}")
assert abs(left - right) < TOL

Compact support: left: -0.138889, right: -0.138889
Integration by parts: left: 0.236111, right: 0.236111
Integration by parts: left: -0.363333, right: -0.363333

The Laplace-Beltrami operator can be expressed by the exterior (co-)derivative

\[\begin{align*} \Delta \alpha = -(d\delta +\delta d)\alpha,\qquad \alpha\in \Lambda^k(\Omega). \end{align*}\]

[10]:

mf = dg.RiemannianManifold(Id(3))

f = dg.ScalarField(20 * x * y * (1 - y) * (1 - z), dim=3)
alpha = dg.OneForm(CF((0.3 * x * y, z**2, -0.1 * x)))

laplace_beltrami = -(dg.d(mf.delta(f)) + mf.delta(dg.d(f)))
expected = mf.CovLaplace(f)
print(f"Delta f: {l2_error(laplace_beltrami, expected, mesh3):.6f}")
assert l2_error(laplace_beltrami, expected, mesh3) < TOL

laplace_beltrami = -(dg.d(mf.delta(alpha)) + mf.delta(dg.d(alpha)))
expected = mf.CovLaplace(alpha)
print(f"Delta alpha: {l2_error(laplace_beltrami, expected, mesh3):.6f}")
assert l2_error(laplace_beltrami, expected, mesh3) < TOL

Delta f: 0.000000
Delta alpha: 0.000000

3D with general Riemannian manifold

Nest, we consider a general Riemannian manifold. Then the Hodge star and exterior co-derivative read in coordinates

\[\begin{align*} . \end{align*}\]

Adapting the volume forms, \(\omega=\sqrt{\det g}\,dx\wedge dy\wedge dz\), \(\omega_{\partial}=\sqrt{\|\mathrm{cof}(g)\|_F}\,dt_1\wedge dt_2\), the \(g\) inner product \(\langle\cdot,\cdot\rangle= g(\cdot,\cdot)\) extended from vector fields to general tensor fields, and the \(g\)-normal vector \(n=n_g\) the above identities hold analogously.

[11]:

import ngsolve


mf = dg.RiemannianManifold(dg.WarpedProduct().metric)

################ covariant div and curl ################
X = dg.VectorField(CF((0.3 * y * z - x * y, x**2 * z + 0.34 * y**3, -x * y * z)))

div_X = mf.star(dg.d(mf.star(mf.Lower(X))))
div_X_expected = mf.CovDiv(X)
print(f"||div X - expected||_L2: {l2_norm(div_X - div_X_expected, mesh3):.6f}")
assert l2_norm(div_X - div_X_expected, mesh3) < TOL

curl_X = mf.Raise(mf.star(dg.d(mf.Lower(X))))
curl_X_expected = mf.CovCurl(X)
print(f"||curl X - expected||_L2: {l2_norm(curl_X - curl_X_expected, mesh3):.6f}")
assert l2_norm(curl_X - curl_X_expected, mesh3) < TOL

################ exterior co-derivative and integration-by-parts ################
f = dg.ScalarField(200 * x * (1 - x) * y * (1 - y) * z * (1 - z), dim=3)
alpha = dg.OneForm(CF((0.3 * x * y, z**2, -0.1 * x)))

left = Integrate(
    mf.InnerProduct(alpha, dg.d(f)) * mf.VolumeForm(VOL) * ngsolve.dx(bonus_intorder=5),
    mesh3,
)
right = Integrate(
    mf.InnerProduct(f, mf.delta(alpha))
    * mf.VolumeForm(VOL)
    * ngsolve.dx(bonus_intorder=5),
    mesh3,
)
print(f"Compact support: left: {left:.6f}, right: {right:.6f}")
assert abs(left - right) < TOL

f = dg.ScalarField(20 * x * y * (1 - y) * (1 - z), dim=3)
alpha = dg.OneForm(CF((0.3 * x * y, z**2, -0.1 * x)))

left = Integrate(
    mf.InnerProduct(alpha, dg.d(f)) * mf.VolumeForm(VOL) * ngsolve.dx(bonus_intorder=5),
    mesh3,
)
right = Integrate(
    mf.InnerProduct(f, mf.delta(alpha))
    * mf.VolumeForm(VOL)
    * ngsolve.dx(bonus_intorder=5),
    mesh3,
) + Integrate(
    f
    * mf.Contraction(mf.normal, alpha)
    * mf.VolumeForm(BND)
    * ngsolve.dx(element_boundary=True, bonus_intorder=5),
    mesh3,
)
print(f"Integration by parts: left: {left:.6f}, right: {right:.6f}")
assert abs(left - right) < TOL

alpha = dg.OneForm(CF((0.3 * x * y, z**2, -0.1 * x)))
beta = dg.TwoForm(
    CF(
        (
            0,
            y * z - 0.1 * x**2,
            0,
            -(y * z - 0.1 * x**2),
            0,
            x * z,
            0,
            -x * z,
            0,
        ),
        dims=(3, 3),
    )
)
left = Integrate(
    0.5
    * mf.InnerProduct(beta, dg.d(alpha))
    * mf.VolumeForm(VOL)
    * ngsolve.dx(bonus_intorder=5),
    mesh3,
)
right = Integrate(
    mf.InnerProduct(alpha, mf.delta(beta))
    * mf.VolumeForm(VOL)
    * ngsolve.dx(bonus_intorder=5),
    mesh3,
) + Integrate(
    mf.InnerProduct(alpha, mf.Contraction(mf.normal, beta), BND)
    * mf.VolumeForm(BND)
    * ngsolve.dx(element_boundary=True, bonus_intorder=5),
    mesh3,
)
print(f"Integration by parts: left: {left:.6f}, right: {right:.6f}")
assert abs(left - right) < TOL

||div X - expected||_L2: 0.000000
||curl X - expected||_L2: 0.000000
Compact support: left: 0.138889, right: 0.138889
Integration by parts: left: 0.479948, right: 0.479948
Integration by parts: left: -0.341228, right: -0.341228

2D Euclidean case

We also present some identities in the two-dimensional Euclidean.

[12]:

mesh2 = Mesh(unit_square.GenerateMesh(maxh=2))

dx = dg.OneForm(CF((1, 0)))
dy = dg.OneForm(CF((0, 1)))

dx_dy = dg.Wedge(dx, dy)

dx_dy2 = dg.TwoForm(CF((0, 1, -1, 0), dims=(2, 2)))

print(f"L2 error dx_dy - dx_dy2 = {l2_error(dx_dy, dx_dy2, mesh2):.6f}")
assert l2_error(dx_dy, dx_dy2, mesh2) < TOL

L2 error dx_dy - dx_dy2 = 0.000000

[13]:

dy_dx = dg.Wedge(dy, dx)
print(f"L2 error dy^dx + dx^dy = {l2_norm(dy_dx + dx_dy, mesh2):.6f}")
assert l2_norm(dy_dx + dx_dy, mesh2) < TOL

L2 error dy^dx + dx^dy = 0.000000

[14]:

mf = dg.RiemannianManifold(Id(2))
one = dg.ScalarField(CF(1), dim=2)
star_one = mf.star(one)
star_dx = mf.star(dx)
star_dy = mf.star(dy)
star_dx_dy = mf.star(dx_dy)
print(f"L2 error star(1) - dx^dy = {l2_error(star_one, dx_dy, mesh2):.6f}")
assert l2_error(star_one, dx_dy, mesh2) < TOL
print(f"L2 error star(dx) - dy = {l2_error(star_dx, dy, mesh2):.6f}")
assert l2_error(star_dx, dy, mesh2) < TOL
print(f"L2 error star(dy) - -dx = {l2_error(star_dy, -dx, mesh2):.6f}")
assert l2_error(star_dy, -dx, mesh2) < TOL
print(f"L2 error star(dx^dy) - 1 = {l2_error(star_dx_dy, one, mesh2):.6f}")
assert l2_error(star_dx_dy, one, mesh2) < TOL

left = dg.Wedge(dx, star_dx)
right = mf.InnerProduct(dx, dx) * dx_dy
print(f"L2 error dx^star(dx) - <dx,dx>*dx^dy = {l2_error(left, right, mesh2):.6f}")
assert l2_error(left, right, mesh2) < TOL
alpha = dg.OneForm(CF((0.3 * x * y, -0.1 * x)))
beta = (y - 0.1 * x**2) * dx_dy

star_star_alpha = mf.star(mf.star(alpha))
expected_star_star_alpha = -alpha
star_star_beta = mf.star(mf.star(beta))
expected_star_star_beta = beta

print(
    f"L2 error star(star(alpha)) + alpha = {l2_error(star_star_alpha, expected_star_star_alpha, mesh2):.6f}"
)
assert l2_error(star_star_alpha, expected_star_star_alpha, mesh2) < TOL
print(
    f"L2 error star(star(beta)) - beta = {l2_error(star_star_beta, expected_star_star_beta, mesh2):.6f}"
)
assert l2_error(star_star_beta, expected_star_star_beta, mesh2) < TOL

L2 error star(1) - dx^dy = 0.000000
L2 error star(dx) - dy = 0.000000
L2 error star(dy) - -dx = 0.000000
L2 error star(dx^dy) - 1 = 0.000000
L2 error dx^star(dx) - <dx,dx>*dx^dy = 0.000000
L2 error star(star(alpha)) + alpha = 0.000000
L2 error star(star(beta)) - beta = 0.000000

[15]:

mesh2 = Mesh(unit_square.GenerateMesh(maxh=0.25))

# scalar 0-form and 1-form
f = dg.ScalarField(x * y, dim=2)
alpha = dg.OneForm(CF((x**2, y**2)))

df = dg.d(f)
ddf = dg.d(df)
print(f"||d^2 f||_L2: {l2_norm(ddf, mesh2):.6f}")
assert l2_norm(ddf, mesh2) < TOL

# Leibniz rule d(f alpha) = df^alpha + f d alpha
left = dg.d(dg.Wedge(f, alpha))
right = dg.Wedge(df, alpha) + dg.Wedge(f, dg.d(alpha))
print(f"Leibniz residual L2: {l2_error(left, right, mesh2):.6f}")
assert l2_error(left, right, mesh2) < TOL

# Alternation of wedge of 1-forms
beta = dg.OneForm(CF((y, x)))
ab = dg.Wedge(alpha, beta)
ba = dg.Wedge(beta, alpha)
print(f"||alpha^beta + beta^alpha||_L2: {l2_norm(ab + ba, mesh2):.6f}")
assert l2_norm(ab + ba, mesh2) < TOL

||d^2 f||_L2: 0.000000
Leibniz residual L2: 0.000000
||alpha^beta + beta^alpha||_L2: 0.000000

[ ]: