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.Einsum("i,j->ij", dx, dy), p=1, q=1, dim=3)
B11 = dg.DoubleForm(dg.Einsum("i,j->ij", alpha1, beta1), p=1, q=1, dim=3)
C11 = dg.DoubleForm(dg.Einsum("i,j->ij", beta1, gamma1), p=1, q=1, dim=3)

A12 = dg.DoubleForm(dg.Einsum("i,jk->ijk", gamma1, alpha2), p=1, q=2, dim=3)
B12 = dg.DoubleForm(dg.Einsum("i,jk->ijk", beta1, beta2), p=1, q=2, dim=3)
C21 = dg.DoubleForm(dg.Einsum("ij,k->ijk", alpha2, dx), p=2, q=1, dim=3)

A22 = dg.Wedge(B11, C11)
B22 = dg.DoubleForm(dg.Einsum("ij,kl->ijkl", alpha2, beta2), p=2, q=2, dim=3)

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

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

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. I.e., the trace contracts the first slot of the left 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\R\) for \(\varphi=\alpha\otimes\beta\), \(\alpha,\beta\in\Lambda^p(\Omega)\) by

\[\begin{align*} \langle A\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

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\gamma) = (\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

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}\varphi^l\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

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