Double forms

In this notebook, we show how to work with \((p,q)\)-double forms by demonstrating several basic identities.

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


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))

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)))

f = dg.ScalarField(20 * x * y * (1 - y) * (1 - z), dim=3)

alpha1 = dg.OneForm(CF((0.3 * x * y, z**2, -0.1 * x)))
beta1 = dg.OneForm(CF((0.3 * z * y, x * z**2, y**2)))
gamma1 = dg.OneForm(CF((0.3 * y * z - x * y, x**2 * z + 0.34 * y**3, -x * y * z)))

alpha2 = dg.TwoForm(
    CF(
        (
            0,
            -0.1 * z**2 - x * y,
            x * y,
            -(-0.1 * z**2 - x * y),
            0,
            0.23 * x * y * z,
            -x * y,
            -0.23 * x * y * z,
            0,
        ),
        dims=(3, 3),
    )
)
beta2 = 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),
    )
)

A11 = dg.DoubleForm(dg.TensorProduct(dx, dy), p=1, q=1, dim=3)
B11 = dg.DoubleForm(dg.TensorProduct(alpha1, beta1), p=1, q=1, dim=3)
C11 = dg.DoubleForm(dg.TensorProduct(beta1, gamma1), p=1, q=1, dim=3)

A12 = dg.DoubleForm(dg.TensorProduct(gamma1, alpha2), p=1, q=2, dim=3)
B12 = dg.DoubleForm(dg.TensorProduct(beta1, beta2), p=1, q=2, dim=3)
C21 = dg.DoubleForm(dg.TensorProduct(alpha2, dx), p=2, q=1, dim=3)

A22 = dg.Wedge(B11, C11)
B22 = dg.DoubleForm(dg.TensorProduct(alpha2, beta2), p=2, q=2, dim=3)

A33 = dg.Wedge(A22, B11)
B33 = dg.Wedge(B22, C11)

Inner products

Define with \(\langle \cdot,\cdot\rangle\) the canonical inner product of \(k\)-forms and double-forms. Further, denote by \(g(\cdot,\cdot)\) the canonical \(g\) inner product of tensor fields. There holds

\[\begin{align*} \langle \varphi,\Psi\rangle=\frac{1}{p!q!}g(\varphi,\Psi),\qquad \varphi,\Psi\in \Lambda^{p,q}(\Omega). \end{align*}\]

# mf = dg.RiemannianManifold(Id(3))
mf = dg.RiemannianManifold(dg.WarpedProduct().metric)

left = mf.InnerProduct(B11, C11)
right = mf.InnerProduct(B11, C11, forms=True)

print(f"L2 error inner product (1,1)-doubleforms: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

left = (
    1
    / (math.factorial(A12.degree_left) * math.factorial(A12.degree_right))
    * mf.InnerProduct(A12, B12)
)
right = mf.InnerProduct(A12, B12, forms=True)
print(f"L2 error inner product (1,2)-doubleforms: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

L2 error inner product (1,1)-doubleforms: 0.000000
L2 error inner product (1,2)-doubleforms: 0.000000

Trace

The trace of a double-form is defined for an orthonormal basis \(\{e_i\}_{i=1}^N\) by

\[\begin{align*} \mathrm{Tr}:\Lambda^{p,q}(\Omega)\to\Lambda^{p-1,q-1}:\quad \mathrm{Tr}(\alpha\otimes\beta)=\sum_{i=1}^N(\iota_{e_i}\alpha)\otimes(\iota_{e_i}\beta), \end{align*}\]

where \(\iota_X\) is the contraction operator for \(k\)-forms. That is, the trace contracts the first slot of the left part with the first slot of the right part of the double form.

We also define the mapping \(\langle \cdot\rangle:\Lambda^{p,p}(\Omega)\to\mathbb{R}\) for \(\varphi=\alpha\otimes\beta\), \(\alpha,\beta\in\Lambda^p(\Omega)\) by

\[\begin{align*} \langle \varphi\rangle=\langle \alpha,\beta\rangle. \end{align*}\]

Define \(\mathrm{Tr}^l:= \mathrm{Tr}\circ\cdots\circ\mathrm{Tr}\) with \(\mathrm{Tr}^0=\mathrm{id}\). There holds for all \(\varphi\in\Lambda^{p,p}(\Omega)\)

\[\begin{align*} \langle \varphi\rangle = \frac{1}{p!}\mathrm{Tr}^p\varphi,\qquad\varphi\in\Lambda^{p,p}(\Omega). \end{align*}\]

left = dg.slot_inner_product(A11, mf)
right = mf.Trace(A11)
print(f"L2 error trace (1,1)-doubleforms: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

left = dg.slot_inner_product(B22, mf, forms=True)
right = 1 / (math.factorial(B22.degree_left)) * mf.Trace(B22, l=2)
print(f"L2 error trace (2,2)-doubleforms: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

L2 error trace (1,1)-doubleforms: 0.000000
L2 error trace (2,2)-doubleforms: 0.000000

Wedge product and Hodge star

Analogously to \(k\)-forms we can define a wedge operator and a Hodge star operator for double-forms, defined as

\[\begin{split}\begin{align*} &\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{\circ\mkern-7mu{\scriptscriptstyle\wedge}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\wedge}\mkern2mu}}:\Lambda^{k,l}(\Omega)\times \Lambda^{p,q}(\Omega)\to\Lambda^{k+p,l+q}(\Omega):\quad&&(\alpha\otimes\beta)\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{\circ\mkern-7mu{\scriptscriptstyle\wedge}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\wedge}\mkern2mu}} (\gamma\otimes\delta) = (\alpha\wedge\gamma)\otimes (\beta\wedge\delta),\\ &\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-9mu{\scriptstyle\star}\mkern4mu}{{\scriptstyle\bigcirc}\mkern-9mu{\scriptstyle\star}\mkern4mu}{\circ\mkern-7mu{\scriptscriptstyle\star}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\star}\mkern2mu}}:\Lambda^{p,q}(\Omega)\to\Lambda^{N-p,N-q}(\Omega):\quad&& \mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-9mu{\scriptstyle\star}\mkern4mu}{{\scriptstyle\bigcirc}\mkern-9mu{\scriptstyle\star}\mkern4mu}{\circ\mkern-7mu{\scriptscriptstyle\star}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\star}\mkern2mu}}(\alpha\otimes\beta) = \star\alpha\otimes\star\beta. \end{align*}\end{split}\]

For example, there holds the relation

\[\begin{align*} \langle \varphi,\Psi\rangle = \mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-9mu{\scriptstyle\star}\mkern4mu}{{\scriptstyle\bigcirc}\mkern-9mu{\scriptstyle\star}\mkern4mu}{\circ\mkern-7mu{\scriptscriptstyle\star}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\star}\mkern2mu}}^{-1}(\varphi\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{\circ\mkern-7mu{\scriptscriptstyle\wedge}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\wedge}\mkern2mu}}\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-9mu{\scriptstyle\star}\mkern4mu}{{\scriptstyle\bigcirc}\mkern-9mu{\scriptstyle\star}\mkern4mu}{\circ\mkern-7mu{\scriptscriptstyle\star}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\star}\mkern2mu}}\Psi),\qquad \varphi,\Psi\in\Lambda^{p,q}(\Omega). \end{align*}\]

left = mf.InnerProduct(C11, B11, forms=True)
right = mf.inv_star(dg.Wedge(C11, mf.star(B11)))
print(f"L2 (1,1)-doubleforms: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

left = mf.InnerProduct(A12, B12, forms=True)
right = mf.inv_star(dg.Wedge(A12, mf.star(B12)))
print(f"L2 (1,2)-doubleforms: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL


left = mf.InnerProduct(A22, B22, forms=True)
right = mf.inv_star(dg.Wedge(A22, mf.star(B22)))
print(f"L2 (2,2)-doubleforms: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

L2 (1,1)-doubleforms: 0.000000
L2 (1,2)-doubleforms: 0.000000
L2 (2,2)-doubleforms: 0.000000

There holds for \(\varphi\in \Lambda^{p,q}(\Omega)\) and \(\Psi\in \Lambda^{p+1,q+1}(\Omega)\)

\[\begin{align*} \langle g\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{\circ\mkern-7mu{\scriptscriptstyle\wedge}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\wedge}\mkern2mu}}\varphi,\Psi\rangle= \langle \varphi,\mathrm{Tr}\Psi\rangle, \end{align*}\]

i.e., \(g\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{\circ\mkern-7mu{\scriptscriptstyle\wedge}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\wedge}\mkern2mu}}\) is adjoint to \(\mathrm{Tr}\).

left = mf.InnerProduct(dg.Wedge(mf.G, B11), B22, forms=True)
right = mf.InnerProduct(B11, mf.Trace(B22), forms=True)
print(
    f"L2 error wedge-trace dual (1,1)-(2,2)-doubleforms: {l2_error(left, right, mesh3):.6f}"
)
assert l2_error(left, right, mesh3) < TOL

A23 = dg.Wedge(A12, B11)
left = mf.InnerProduct(dg.Wedge(mf.G, B12), A23, forms=True)
right = mf.InnerProduct(B12, mf.Trace(A23), forms=True)
print(
    f"L2 error wedge-trace dual (1,2)-(2,3)-doubleforms: {l2_error(left, right, mesh3):.6f}"
)
assert l2_error(left, right, mesh3) < TOL

L2 error wedge-trace dual (1,1)-(2,2)-doubleforms: 0.000000
L2 error wedge-trace dual (1,2)-(2,3)-doubleforms: 0.000000

Metric powers and contraction identities

Let \(g^l=g\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{\circ\mkern-7mu{\scriptscriptstyle\wedge}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\wedge}\mkern2mu}}\cdots\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{\circ\mkern-7mu{\scriptscriptstyle\wedge}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\wedge}\mkern2mu}} g\) with \(g^0:=1\). There holds for all \(\varphi\in \Lambda^{p,q}(\Omega)\) and \(0\le l\le\min\{N-p,N-q\}\)

\[\begin{align*} \mathrm{Tr}^l\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-9mu{\scriptstyle\star}\mkern4mu}{{\scriptstyle\bigcirc}\mkern-9mu{\scriptstyle\star}\mkern4mu}{\circ\mkern-7mu{\scriptscriptstyle\star}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\star}\mkern2mu}}\varphi = \mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-9mu{\scriptstyle\star}\mkern4mu}{{\scriptstyle\bigcirc}\mkern-9mu{\scriptstyle\star}\mkern4mu}{\circ\mkern-7mu{\scriptscriptstyle\star}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\star}\mkern2mu}}(\varphi\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{\circ\mkern-7mu{\scriptscriptstyle\wedge}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\wedge}\mkern2mu}} g^l). \end{align*}\]

G2 = dg.Wedge(mf.G, mf.G)
G3 = dg.Wedge(G2, mf.G)

left = mf.Trace(mf.star(B11), l=1)
right = mf.star(dg.Wedge(B11, mf.G))
print(f"L2 error (1,1)-doubleforms l=1: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

left = mf.Trace(mf.star(B11), l=2)
right = mf.star(dg.Wedge(B11, G2))
print(f"L2 error (1,1)-doubleforms l=2: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

left = mf.Trace(mf.star(B22), l=1)
right = mf.star(dg.Wedge(B22, mf.G))
print(f"L2 error (2,2)-doubleforms l=1: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

left = mf.Trace(mf.star(C21), l=1)
right = mf.star(dg.Wedge(C21, mf.G))
print(f"L2 error (2,1)-doubleforms l=1: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

left = mf.Trace(mf.star(B12), l=1)
right = mf.star(dg.Wedge(B12, mf.G))
print(f"L2 error (1,2)-doubleforms l=1: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

L2 error (1,1)-doubleforms l=1: 0.000000
L2 error (1,1)-doubleforms l=2: 0.000000
L2 error (2,2)-doubleforms l=1: 0.000000
L2 error (2,1)-doubleforms l=1: 0.000000
L2 error (1,2)-doubleforms l=1: 0.000000

For all \(\varphi\in \Lambda^{p,p}(\Omega)\) and \(0\le l \le N-p\) there holds

\[\begin{align*} \langle \varphi\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{\circ\mkern-7mu{\scriptscriptstyle\wedge}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\wedge}\mkern2mu}} g^l\rangle = \frac{(N-p)!}{(N-p-l)!}\langle \varphi\rangle. \end{align*}\]

def right_term(phi, l, N=3):
    p = phi.degree_left
    return (
        math.factorial(N - p)
        / math.factorial(N - p - l)
        * dg.slot_inner_product(phi, mf)
    )


left = dg.slot_inner_product(dg.Wedge(B11, mf.G), mf)
right = right_term(B11, 1)
print(f"L2 error (1,1)-doubleforms l=1: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

left = dg.slot_inner_product(dg.Wedge(B11, dg.Wedge(mf.G, mf.G)), mf)
right = right_term(B11, 2)
print(f"L2 error (1,1)-doubleforms l=2: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL


left = dg.slot_inner_product(dg.Wedge(A22, mf.G), mf)
right = right_term(A22, 1)
print(f"L2 error (2,2)-doubleforms l=1: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

L2 error (1,1)-doubleforms l=1: 0.000000
L2 error (1,1)-doubleforms l=2: 0.000000
L2 error (2,2)-doubleforms l=1: 0.000000

For all \(\varphi\in\Lambda^{p,p}(\Omega)\) there holds

\[\begin{align*} \mathrm{Tr}(\varphi\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{\circ\mkern-7mu{\scriptscriptstyle\wedge}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\wedge}\mkern2mu}} g) = (N-2p)\varphi+(\mathrm{Tr}\varphi)\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{\circ\mkern-7mu{\scriptscriptstyle\wedge}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\wedge}\mkern2mu}} g. \end{align*}\]

There holds for all \(\varphi\in\Lambda^{p,p}(\Omega)\) and \(l\ge 1\)

\[\begin{align*} \mathrm{Tr}(\varphi\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{\circ\mkern-7mu{\scriptscriptstyle\wedge}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\wedge}\mkern2mu}} g^l)=l(N-2p-l+1)\varphi\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{\circ\mkern-7mu{\scriptscriptstyle\wedge}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\wedge}\mkern2mu}} g^{l-1}+\mathrm{Tr}\varphi\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{\circ\mkern-7mu{\scriptscriptstyle\wedge}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\wedge}\mkern2mu}} g^l. \end{align*}\]

def right_term(phi, l, N=3):
    p = phi.degree_left
    return l * (N - 2 * p - l + 1) * dg.Wedge(
        phi, dg.WedgePower(mf.G, l - 1)
    ) + dg.Wedge(mf.Trace(phi), dg.WedgePower(mf.G, l))


left = mf.Trace(dg.Wedge(B11, mf.G))
right = right_term(B11, 1)
print(f"L2 error (1,1)-doubleforms l=1: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

left = mf.Trace(dg.Wedge(B11, dg.WedgePower(mf.G, 2)))
right = right_term(B11, 2)
print(f"L2 error (1,1)-doubleforms l=2: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

left = mf.Trace(dg.Wedge(A22, mf.G))
right = right_term(A22, 1)
print(f"L2 error (2,2)-doubleforms l=1: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

left = mf.Trace(dg.Wedge(A22, dg.WedgePower(mf.G, 2)))
right = right_term(A22, 2)
print(f"L2 error (2,2)-doubleforms l=2: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

L2 error (1,1)-doubleforms l=1: 0.000000
L2 error (1,1)-doubleforms l=2: 0.000000
L2 error (2,2)-doubleforms l=1: 0.000000
L2 error (2,2)-doubleforms l=2: 0.000000

For all \(\varphi\in\Lambda^{p,p}(\Omega)\) and all \(l\ge 1\) there holds

\[\begin{align*} \mathrm{Tr}^l(\varphi\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{\circ\mkern-7mu{\scriptscriptstyle\wedge}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\wedge}\mkern2mu}} g) = l(N-2p+l-1)\mathrm{Tr}^{l-1}\varphi+(\mathrm{Tr}^l\varphi)\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{\circ\mkern-7mu{\scriptscriptstyle\wedge}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\wedge}\mkern2mu}} g. \end{align*}\]

def right_term(phi, l, N=3):
    p = phi.degree_left
    return l * (N - 2 * p + l - 1) * mf.Trace(phi, l=l - 1) + dg.Wedge(
        mf.Trace(phi, l=l), mf.G
    )


left = mf.Trace(dg.Wedge(B11, mf.G), l=1)
right = right_term(B11, 1)
print(f"L2 error (1,1)-doubleforms l=1: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

left = mf.Trace(dg.Wedge(B11, mf.G), l=2)
right = right_term(B11, 2)
print(f"L2 error (1,1)-doubleforms l=2: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

left = mf.Trace(dg.Wedge(A22, mf.G), l=1)
right = right_term(A22, 1)
print(f"L2 error (2,2)-doubleforms l=1: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

left = mf.Trace(dg.Wedge(A22, mf.G), l=2)
right = right_term(A22, 2)
print(f"L2 error (2,2)-doubleforms l=2: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

L2 error (1,1)-doubleforms l=1: 0.000000
L2 error (1,1)-doubleforms l=2: 0.000000
L2 error (2,2)-doubleforms l=1: 0.000000
L2 error (2,2)-doubleforms l=2: 0.000000

For all \(p\ge l\ge0\) there holds

\[\begin{align*} \frac{1}{l!}\mathrm{Tr}^l\left(\frac{1}{p!}g^p\right) = \binom{N-p+l}{l}\frac{1}{(p-l)!}g^{p-l}. \end{align*}\]

def term_left(p, l):
    if l < 0 or p < l:
        raise ValueError("Invalid values for p and l.")
    return (
        1
        / math.factorial(l)
        * mf.Trace(1 / math.factorial(p) * dg.WedgePower(mf.G, p), l=l)
    )


def term_right(p, l, N=3):
    if l < 0 or p < l:
        raise ValueError("Invalid values for p and l.")
    return (
        math.comb(N - p + l, l) * 1 / math.factorial(p - l) * dg.WedgePower(mf.G, p - l)
    )


for p in range(3):
    for l in range(p + 1):
        left = term_left(p, l)
        right = term_right(p, l)
        print(f"L2 error p={p}, l={l}: {l2_error(left, right, mesh3):.6f}")
        assert l2_error(left, right, mesh3) < TOL

L2 error p=0, l=0: 0.000000
L2 error p=1, l=0: 0.000000
L2 error p=1, l=1: 0.000000
L2 error p=2, l=0: 0.000000
L2 error p=2, l=1: 0.000000
L2 error p=2, l=2: 0.000000

There holds with the volume form \(\omega\in\Lambda^N(\Omega)\)

\[\begin{align*} \frac{1}{N!}g^N = \omega\otimes\omega. \end{align*}\]

omega = mf.VolumeForm(VOL) * dg.Wedge(dx, dg.Wedge(dy, dz))

left = 1 / math.factorial(3) * dg.WedgePower(mf.G, 3)
right = dg.TensorProduct(omega, omega)
print(f"L2 error volume form: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

L2 error volume form: 0.000000

There holds

\[\begin{align*} \frac{1}{(N-1)!}\mathrm{Tr}^{N-1}(\omega\otimes\omega) = g. \end{align*}\]

omega = mf.VolumeForm(VOL) * dg.Wedge(dx, dg.Wedge(dy, dz))
N = 3

left = (
    1
    / math.factorial(N - 1)
    * mf.Trace(dg.DoubleForm(dg.TensorProduct(omega, omega), p=3, q=3, dim=3), l=N - 1)
)
right = mf.G
print(f"L2 error metric from volume form: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

L2 error metric from volume form: 0.000000

There holds for all \(\varphi,\Psi\in\Lambda^{p,p}(\Omega)\)

\[\begin{align*} \langle \varphi\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{\circ\mkern-7mu{\scriptscriptstyle\wedge}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\wedge}\mkern2mu}} g,\Psi\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{\circ\mkern-7mu{\scriptscriptstyle\wedge}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\wedge}\mkern2mu}} g\rangle = (N-2p)\langle \varphi,\Psi\rangle + \langle \mathrm{Tr}\varphi,\mathrm{Tr}\Psi\rangle. \end{align*}\]

N = 3


def left_term(phi, psi):
    return mf.InnerProduct(dg.Wedge(phi, mf.G), dg.Wedge(psi, mf.G), forms=True)


def right_term(phi, psi):
    p = phi.degree_left
    return (N - 2 * p) * mf.InnerProduct(phi, psi, forms=True) + mf.InnerProduct(
        mf.Trace(phi), mf.Trace(psi), forms=True
    )


left = left_term(B11, C11)
right = right_term(B11, C11)
print(f"L2 error (1,1)-doubleforms via formula: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

left = left_term(A22, B22)
right = right_term(A22, B22)
print(f"L2 error (2,2)-doubleforms via formula: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL


left = left_term(A33, B33)
right = right_term(A33, B33)
print(f"L2 error (3,3)-doubleforms via formula: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

L2 error (1,1)-doubleforms via formula: 0.000000
L2 error (2,2)-doubleforms via formula: 0.000000
L2 error (3,3)-doubleforms via formula: 0.000000

Algebraic operators

Define \(\mathcal{S}:\Lambda^{1,1}(\Omega)\to\Lambda^{1,1}(\Omega)\) and \(\mathcal{J}:\Lambda^{1,1}(\Omega)\to\Lambda^{1,1}(\Omega)\) by

\[\begin{split}\begin{align*} &\mathcal{S}\sigma = \sigma^T-g\,\mathrm{Tr}\sigma,\\ &\mathcal{J}\sigma=\sigma^T-\frac{1}{2}g\,\mathrm{Tr}\sigma. \end{align*}\end{split}\]

For all \(\sigma\in\Lambda^{1,1}(\Omega)\) there holds

\[\begin{split}\begin{align*} &\mathcal{S}\sigma = \frac{-1}{(N-2)!}\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-9mu{\scriptstyle\star}\mkern4mu}{{\scriptstyle\bigcirc}\mkern-9mu{\scriptstyle\star}\mkern4mu}{\circ\mkern-7mu{\scriptscriptstyle\star}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\star}\mkern2mu}}^{-1}(g^{N-2}\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{\circ\mkern-7mu{\scriptscriptstyle\wedge}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\wedge}\mkern2mu}}\sigma),\\ &g\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{\circ\mkern-7mu{\scriptscriptstyle\wedge}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\wedge}\mkern2mu}}\mathcal{J}\sigma = \frac{-1}{(N-3)!}\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-9mu{\scriptstyle\star}\mkern4mu}{{\scriptstyle\bigcirc}\mkern-9mu{\scriptstyle\star}\mkern4mu}{\circ\mkern-7mu{\scriptscriptstyle\star}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\star}\mkern2mu}}^{-1}(g^{N-3}\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{\circ\mkern-7mu{\scriptscriptstyle\wedge}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\wedge}\mkern2mu}}\sigma). \end{align*}\end{split}\]

N = 3

left = mf.S(B11)
right = -1 / math.factorial(N - 2) * mf.inv_star(dg.Wedge(mf.G, B11))
print(f"L2 error B operator (1,1)-doubleforms: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

left = dg.Wedge(mf.G, mf.J(B11))
right = -1 / math.factorial(N - 3) * mf.inv_star(B11)
print(f"L2 error J operator (1,1)-doubleforms: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

L2 error B operator (1,1)-doubleforms: 0.000000
L2 error J operator (1,1)-doubleforms: 0.000000

Define the operator \(s:\Lambda^{p,q}(\Omega)\to\Lambda^{p+1,q-1}(\Omega)\) with an orthonormal basis \(\{e_i\}_{i=1}^N\) and its associated dual basis \(\{e^i\}_{i=1}^N\) via

\[\begin{align*} s(\alpha\otimes\beta)=\sum_{i=1}^N(e^i\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{\circ\mkern-7mu{\scriptscriptstyle\wedge}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\wedge}\mkern2mu}}\alpha)\otimes (\iota_{e_i}\beta). \end{align*}\]

The maps \(\mathrm{Tr}\) and \(s\) commute.

left = mf.Trace(mf.s(A22))
right = mf.s(mf.Trace(A22))
print(
    f"L2 error trace-s operator (2,2)-doubleforms: {l2_error(left, right, mesh3):.6f}"
)
assert l2_error(left, right, mesh3) < TOL

left = mf.Trace(mf.s(A33))
right = mf.s(mf.Trace(A33))
print(
    f"L2 error trace-s operator (3,3)-doubleforms: {l2_error(left, right, mesh3):.6f}"
)
assert l2_error(left, right, mesh3) < TOL

L2 error trace-s operator (2,2)-doubleforms: 0.000000
L2 error trace-s operator (3,3)-doubleforms: 0.000000

There holds the Leibniz rule for \(s\) for all \(\varphi\in\Lambda^{p,q}(\Omega)\) and \(\Psi\in\Lambda^{l,k}(\Omega)\)

\[\begin{align*} s(\varphi\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{\circ\mkern-7mu{\scriptscriptstyle\wedge}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\wedge}\mkern2mu}}\Psi)=(s\varphi)\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{\circ\mkern-7mu{\scriptscriptstyle\wedge}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\wedge}\mkern2mu}}\Psi + (-1)^{p+q}\varphi\mathbin{\mathchoice{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{{\scriptstyle\bigcirc}\mkern-10mu{\scriptstyle\wedge}\mkern3mu}{\circ\mkern-7mu{\scriptscriptstyle\wedge}\mkern3mu}{\circ\mkern-6mu{\scriptscriptstyle\wedge}\mkern2mu}}(s\Psi). \end{align*}\]

def left_term(phi, psi):
    return mf.s(dg.Wedge(phi, psi))


def right_term(phi, psi):
    p = phi.degree_left
    q = phi.degree_right
    return dg.Wedge(mf.s(phi), psi) + (-1) ** (p + q) * dg.Wedge(phi, mf.s(psi))


left = left_term(B11, C11)
right = right_term(B11, C11)
print(f"L2 error (1,1)-(1,1)-doubleforms: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

left = left_term(A22, B12)
right = right_term(A22, B12)
print(f"L2 error (2,2)-(1,2)-doubleforms: {l2_error(left, right, mesh3):.6f}")
assert l2_error(left, right, mesh3) < TOL

L2 error (1,1)-(1,1)-doubleforms: 0.000000
L2 error (2,2)-(1,2)-doubleforms: 0.000000