{ "cells": [ { "cell_type": "markdown", "id": "3995b4f7", "metadata": {}, "source": [ "# Double-forms\n", "In this notebook, we show how to work with $(p,q)$-double forms by demonstrating several basic identities.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "ed476010", "metadata": {}, "outputs": [], "source": [ "from ngsolve import *\n", "from netgen.occ import unit_square, unit_cube\n", "import ngsdiffgeo as dg\n", "import math\n", "\n", "\n", "TOL = 1e-8\n", "\n", "\n", "def l2_error(a, b, mesh):\n", " return sqrt(Integrate(InnerProduct(a - b, a - b), mesh))\n", "\n", "\n", "def l2_norm(a, mesh):\n", " return sqrt(Integrate(InnerProduct(a, a), mesh))" ] }, { "cell_type": "code", "execution_count": null, "id": "68c140b1", "metadata": {}, "outputs": [], "source": [ "mesh3 = Mesh(unit_cube.GenerateMesh(maxh=2))\n", "\n", "dx = dg.OneForm(CF((1, 0, 0)))\n", "dy = dg.OneForm(CF((0, 1, 0)))\n", "dz = dg.OneForm(CF((0, 0, 1)))\n", "\n", "f = dg.ScalarField(20 * x * y * (1 - y) * (1 - z), dim=3)\n", "\n", "alpha1 = dg.OneForm(CF((0.3 * x * y, z**2, -0.1 * x)))\n", "beta1 = dg.OneForm(CF((0.3 * z * y, x * z**2, y**2)))\n", "gamma1 = dg.OneForm(CF((0.3 * y * z - x * y, x**2 * z + 0.34 * y**3, -x * y * z)))\n", "\n", "alpha2 = dg.TwoForm(\n", " CF(\n", " (\n", " 0,\n", " -0.1 * z**2 - x * y,\n", " x * y,\n", " -(-0.1 * z**2 - x * y),\n", " 0,\n", " 0.23 * x * y * z,\n", " -x * y,\n", " -0.23 * x * y * z,\n", " 0,\n", " ),\n", " dims=(3, 3),\n", " )\n", ")\n", "beta2 = dg.TwoForm(\n", " CF(\n", " (\n", " 0,\n", " y * z - 0.1 * x**2,\n", " 0,\n", " -(y * z - 0.1 * x**2),\n", " 0,\n", " x * z,\n", " 0,\n", " -x * z,\n", " 0,\n", " ),\n", " dims=(3, 3),\n", " )\n", ")\n", "\n", "A11 = dg.DoubleForm(dg.Einsum(\"i,j->ij\", dx, dy), p=1, q=1, dim=3)\n", "B11 = dg.DoubleForm(dg.Einsum(\"i,j->ij\", alpha1, beta1), p=1, q=1, dim=3)\n", "C11 = dg.DoubleForm(dg.Einsum(\"i,j->ij\", beta1, gamma1), p=1, q=1, dim=3)\n", "\n", "A12 = dg.DoubleForm(dg.Einsum(\"i,jk->ijk\", gamma1, alpha2), p=1, q=2, dim=3)\n", "B12 = dg.DoubleForm(dg.Einsum(\"i,jk->ijk\", beta1, beta2), p=1, q=2, dim=3)\n", "C21 = dg.DoubleForm(dg.Einsum(\"ij,k->ijk\", alpha2, dx), p=2, q=1, dim=3)\n", "\n", "A22 = dg.Wedge(B11, C11)\n", "B22 = dg.DoubleForm(dg.Einsum(\"ij,kl->ijkl\", alpha2, beta2), p=2, q=2, dim=3)\n", "\n", "A33 = dg.Wedge(A22, B11)\n", "B33 = dg.Wedge(B22, C11)\n" ] }, { "cell_type": "markdown", "id": "bd219589", "metadata": {}, "source": [ "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\n", "\n", "$$\n", "\\begin{align*}\n", "\\langle \\varphi,\\Psi\\rangle=\\frac{1}{p!q!}g(\\varphi,\\Psi),\\qquad \\varphi,\\Psi\\in \\Lambda^{p,q}(\\Omega).\n", "\\end{align*}\n", "$$\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "28d677cc", "metadata": {}, "outputs": [], "source": [ "# mf = dg.RiemannianManifold(Id(3))\n", "mf = dg.RiemannianManifold(dg.WarpedProduct().metric)\n", "\n", "left = mf.InnerProduct(B11, C11)\n", "right = mf.InnerProduct(B11, C11, forms=True)\n", "\n", "print(f\"L2 error inner product (1,1)-doubleforms: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL\n", "\n", "left = (\n", " 1\n", " / (math.factorial(A12.degree_left) * math.factorial(A12.degree_right))\n", " * mf.InnerProduct(A12, B12)\n", ")\n", "right = mf.InnerProduct(A12, B12, forms=True)\n", "print(f\"L2 error inner product (1,2)-doubleforms: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL" ] }, { "cell_type": "markdown", "id": "75d35ebe", "metadata": {}, "source": [ "The trace of a double-form is defined for an orthonormal basis $\\{e_i\\}_{i=1}^N$ by\n", "\n", "$$\n", "\\begin{align*}\n", "\\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),\n", "\\end{align*}\n", "$$\n", "\n", "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.\n", "\n", "We also define the mapping\n", "$\\langle \\cdot\\rangle:\\Lambda^{p,p}(\\Omega)\\to\\R$ for $\\varphi=\\alpha\\otimes\\beta$, $\\alpha,\\beta\\in\\Lambda^p(\\Omega)$ by \n", "\n", "$$\n", "\\begin{align*}\n", "\\langle A\\rangle=\\langle \\alpha,\\beta\\rangle.\n", "\\end{align*}\n", "$$\n", "\n", "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)$\n", "\n", "$$\n", "\\begin{align*}\n", "\\langle \\varphi\\rangle = \\frac{1}{p!}\\mathrm{Tr}^p\\varphi,\\qquad\\varphi\\in\\Lambda^{p,p}(\\Omega).\n", "\\end{align*}\n", "$$\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "0e8cae73", "metadata": {}, "outputs": [], "source": [ "left = dg.slot_inner_product(A11, mf)\n", "right = mf.Trace(A11)\n", "print(f\"L2 error trace (1,1)-doubleforms: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL\n", "\n", "left = dg.slot_inner_product(B22, mf, forms=True)\n", "right = 1 / (math.factorial(B22.degree_left)) * mf.Trace(B22, l=2)\n", "print(f\"L2 error trace (2,2)-doubleforms: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL" ] }, { "cell_type": "markdown", "id": "ce7eca05", "metadata": {}, "source": [ "Analogously to $k$-forms we can define a wedge operator and a Hodge star operator for double-forms, defined as\n", "\n", "$$\n", "\\begin{align*}\n", "&\\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),\\\\\n", "&\\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.\n", "\\end{align*}\n", "$$\n", "\n", "For example, there holds the relation\n", "\n", "$$\n", "\\begin{align*}\n", "\\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).\n", "\\end{align*}\n", "$$\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "d5ca8d7a", "metadata": {}, "outputs": [], "source": [ "left = mf.InnerProduct(C11, B11, forms=True)\n", "right = mf.inv_star(dg.Wedge(C11, mf.star(B11)))\n", "print(f\"L2 (1,1)-doubleforms: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL\n", "\n", "left = mf.InnerProduct(A12, B12, forms=True)\n", "right = mf.inv_star(dg.Wedge(A12, mf.star(B12)))\n", "print(f\"L2 (1,2)-doubleforms: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL\n", "\n", "\n", "left = mf.InnerProduct(A22, B22, forms=True)\n", "right = mf.inv_star(dg.Wedge(A22, mf.star(B22)))\n", "print(f\"L2 (2,2)-doubleforms: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL" ] }, { "cell_type": "markdown", "id": "8fce0af4", "metadata": {}, "source": [ "There holds for $\\varphi\\in \\Lambda^{p,q}(\\Omega)$ and $\\Psi\\in \\Lambda^{p+1,q+1}(\\Omega)$\n", "\n", "$$\n", "\\begin{align*}\n", "\\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,\n", "\\end{align*}\n", "$$\n", "\n", "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}$.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "abe7e00c", "metadata": {}, "outputs": [], "source": [ "left = mf.InnerProduct(dg.Wedge(mf.G, B11), B22, forms=True)\n", "right = mf.InnerProduct(B11, mf.Trace(B22), forms=True)\n", "print(\n", " f\"L2 error wedge-trace dual (1,1)-(2,2)-doubleforms: {l2_error(left, right, mesh3):.6f}\"\n", ")\n", "assert l2_error(left, right, mesh3) < TOL\n", "\n", "A23 = dg.Wedge(A12, B11)\n", "left = mf.InnerProduct(dg.Wedge(mf.G, B12), A23, forms=True)\n", "right = mf.InnerProduct(B12, mf.Trace(A23), forms=True)\n", "print(\n", " f\"L2 error wedge-trace dual (1,2)-(2,3)-doubleforms: {l2_error(left, right, mesh3):.6f}\"\n", ")\n", "assert l2_error(left, right, mesh3) < TOL" ] }, { "cell_type": "markdown", "id": "949833c0", "metadata": {}, "source": [ "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\\}$\n", "\n", "$$\n", "\\begin{align*}\n", "\\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).\n", "\\end{align*}\n", "$$\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "0e56bcd2", "metadata": {}, "outputs": [], "source": [ "G2 = dg.Wedge(mf.G, mf.G)\n", "G3 = dg.Wedge(G2, mf.G)\n", "\n", "left = mf.Trace(mf.star(B11), l=1)\n", "right = mf.star(dg.Wedge(B11, mf.G))\n", "print(f\"L2 error (1,1)-doubleforms l=1: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL\n", "\n", "left = mf.Trace(mf.star(B11), l=2)\n", "right = mf.star(dg.Wedge(B11, G2))\n", "print(f\"L2 error (1,1)-doubleforms l=2: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL\n", "\n", "left = mf.Trace(mf.star(B22), l=1)\n", "right = mf.star(dg.Wedge(B22, mf.G))\n", "print(f\"L2 error (2,2)-doubleforms l=1: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL\n", "\n", "left = mf.Trace(mf.star(C21), l=1)\n", "right = mf.star(dg.Wedge(C21, mf.G))\n", "print(f\"L2 error (2,1)-doubleforms l=1: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL\n", "\n", "left = mf.Trace(mf.star(B12), l=1)\n", "right = mf.star(dg.Wedge(B12, mf.G))\n", "print(f\"L2 error (1,2)-doubleforms l=1: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL" ] }, { "cell_type": "markdown", "id": "180aadcd", "metadata": {}, "source": [ "For all $\\varphi\\in \\Lambda^{p,p}(\\Omega)$ and $0\\le l \\le N-p$ there holds\n", "\n", "$$\n", "\\begin{align*}\n", "\\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.\n", "\\end{align*}\n", "$$\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "eccb02bb", "metadata": {}, "outputs": [], "source": [ "def right_term(phi, l, N=3):\n", " p = phi.degree_left\n", " return (\n", " math.factorial(N - p)\n", " / math.factorial(N - p - l)\n", " * dg.slot_inner_product(phi, mf)\n", " )\n", "\n", "\n", "left = dg.slot_inner_product(dg.Wedge(B11, mf.G), mf)\n", "right = right_term(B11, 1)\n", "print(f\"L2 error (1,1)-doubleforms l=1: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL\n", "\n", "left = dg.slot_inner_product(dg.Wedge(B11, dg.Wedge(mf.G, mf.G)), mf)\n", "right = right_term(B11, 2)\n", "print(f\"L2 error (1,1)-doubleforms l=2: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL\n", "\n", "\n", "left = dg.slot_inner_product(dg.Wedge(A22, mf.G), mf)\n", "right = right_term(A22, 1)\n", "print(f\"L2 error (2,2)-doubleforms l=1: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL" ] }, { "cell_type": "markdown", "id": "f912620d", "metadata": {}, "source": [ "For all $\\varphi\\in\\Lambda^{p,p}(\\Omega)$ there holds\n", "\n", "$$\n", "\\begin{align*}\n", "\\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.\n", "\\end{align*}\n", "$$\n", "\n", "There holds for all $\\varphi\\in\\Lambda^{p,p}(\\Omega)$ and $l\\ge 1$\n", "\n", "$$\n", "\\begin{align*}\n", "\\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.\n", "\\end{align*}\n", "$$\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "8e1c7cae", "metadata": {}, "outputs": [], "source": [ "def right_term(phi, l, N=3):\n", " p = phi.degree_left\n", " return l * (N - 2 * p - l + 1) * dg.Wedge(\n", " phi, dg.WedgePower(mf.G, l - 1)\n", " ) + dg.Wedge(mf.Trace(phi), dg.WedgePower(mf.G, l))\n", "\n", "\n", "left = mf.Trace(dg.Wedge(B11, mf.G))\n", "right = right_term(B11, 1)\n", "print(f\"L2 error (1,1)-doubleforms l=1: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL\n", "\n", "left = mf.Trace(dg.Wedge(B11, dg.WedgePower(mf.G, 2)))\n", "right = right_term(B11, 2)\n", "print(f\"L2 error (1,1)-doubleforms l=2: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL\n", "\n", "left = mf.Trace(dg.Wedge(A22, mf.G))\n", "right = right_term(A22, 1)\n", "print(f\"L2 error (2,2)-doubleforms l=1: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL\n", "\n", "left = mf.Trace(dg.Wedge(A22, dg.WedgePower(mf.G, 2)))\n", "right = right_term(A22, 2)\n", "print(f\"L2 error (2,2)-doubleforms l=2: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL" ] }, { "cell_type": "markdown", "id": "5db78c1a", "metadata": {}, "source": [ "For all $\\varphi\\in\\Lambda^{p,p}(\\Omega)$ and all $l\\ge 1$ there holds\n", "\n", "$$\n", "\\begin{align*}\n", "\\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.\n", "\\end{align*}\n", "$$\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "e6be03b3", "metadata": {}, "outputs": [], "source": [ "def right_term(phi, l, N=3):\n", " p = phi.degree_left\n", " return l * (N - 2 * p + l - 1) * mf.Trace(phi, l=l - 1) + dg.Wedge(\n", " mf.Trace(phi, l=l), mf.G\n", " )\n", "\n", "\n", "left = mf.Trace(dg.Wedge(B11, mf.G), l=1)\n", "right = right_term(B11, 1)\n", "print(f\"L2 error (1,1)-doubleforms l=1: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL\n", "\n", "left = mf.Trace(dg.Wedge(B11, mf.G), l=2)\n", "right = right_term(B11, 2)\n", "print(f\"L2 error (1,1)-doubleforms l=2: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL\n", "\n", "left = mf.Trace(dg.Wedge(A22, mf.G), l=1)\n", "right = right_term(A22, 1)\n", "print(f\"L2 error (2,2)-doubleforms l=1: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL\n", "\n", "left = mf.Trace(dg.Wedge(A22, mf.G), l=2)\n", "right = right_term(A22, 2)\n", "print(f\"L2 error (2,2)-doubleforms l=2: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL" ] }, { "cell_type": "markdown", "id": "e7003bb3", "metadata": {}, "source": [ "For all $p\\ge l\\ge0$ there holds\n", "\n", "$$\n", "\\begin{align*}\n", "\\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}.\n", "\\end{align*}\n", "$$\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "62ef91e8", "metadata": {}, "outputs": [], "source": [ "def term_left(p, l):\n", " if l < 0 or p < l:\n", " raise ValueError(\"Invalid values for p and l.\")\n", " return (\n", " 1\n", " / math.factorial(l)\n", " * mf.Trace(1 / math.factorial(p) * dg.WedgePower(mf.G, p), l=l)\n", " )\n", "\n", "\n", "def term_right(p, l, N=3):\n", " if l < 0 or p < l:\n", " raise ValueError(\"Invalid values for p and l.\")\n", " return (\n", " math.comb(N - p + l, l) * 1 / math.factorial(p - l) * dg.WedgePower(mf.G, p - l)\n", " )\n", "\n", "\n", "for p in range(3):\n", " for l in range(p + 1):\n", " left = term_left(p, l)\n", " right = term_right(p, l)\n", " print(f\"L2 error p={p}, l={l}: {l2_error(left, right, mesh3):.6f}\")\n", " assert l2_error(left, right, mesh3) < TOL" ] }, { "cell_type": "markdown", "id": "c413f12c", "metadata": {}, "source": [ "There holds with the volume form $\\omega\\in\\Lambda^N(\\Omega)$\n", "\n", "$$\n", "\\begin{align*}\n", "\\frac{1}{N!}g^N = \\omega\\otimes\\omega.\n", "\\end{align*}\n", "$$\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "38c49723", "metadata": {}, "outputs": [], "source": [ "omega = mf.VolumeForm(VOL) * dg.Wedge(dx, dg.Wedge(dy, dz))\n", "\n", "left = 1 / math.factorial(3) * dg.WedgePower(mf.G, 3)\n", "right = dg.TensorProduct(omega, omega)\n", "print(f\"L2 error volume form: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL" ] }, { "cell_type": "markdown", "id": "f91775d5", "metadata": {}, "source": [ "There holds\n", "\n", "$$\n", "\\begin{align*}\n", "\\frac{1}{(N-1)!}\\mathrm{Tr}^{N-1}(\\omega\\otimes\\omega) = g.\n", "\\end{align*}\n", "$$\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "d38f51ab", "metadata": {}, "outputs": [], "source": [ "omega = mf.VolumeForm(VOL) * dg.Wedge(dx, dg.Wedge(dy, dz))\n", "N = 3\n", "\n", "left = (\n", " 1\n", " / math.factorial(N - 1)\n", " * mf.Trace(dg.DoubleForm(dg.TensorProduct(omega, omega), p=3, q=3, dim=3), l=N - 1)\n", ")\n", "right = mf.G\n", "print(f\"L2 error metric from volume form: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL" ] }, { "cell_type": "markdown", "id": "c8e0cdf5", "metadata": {}, "source": [ "There holds for all $\\varphi,\\Psi\\in\\Lambda^{p,p}(\\Omega)$\n", "\n", "$$\n", "\\begin{align*}\n", "\\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.\n", "\\end{align*}\n", "$$\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "45492b2c", "metadata": {}, "outputs": [], "source": [ "N = 3\n", "\n", "\n", "def left_term(phi, psi):\n", " return mf.InnerProduct(dg.Wedge(phi, mf.G), dg.Wedge(psi, mf.G), forms=True)\n", "\n", "\n", "def right_term(phi, psi):\n", " p = phi.degree_left\n", " return (N - 2 * p) * mf.InnerProduct(phi, psi, forms=True) + mf.InnerProduct(\n", " mf.Trace(phi), mf.Trace(psi), forms=True\n", " )\n", "\n", "\n", "left = left_term(B11, C11)\n", "right = right_term(B11, C11)\n", "print(f\"L2 error (1,1)-doubleforms via formula: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL\n", "\n", "left = left_term(A22, B22)\n", "right = right_term(A22, B22)\n", "print(f\"L2 error (2,2)-doubleforms via formula: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL\n", "\n", "\n", "left = left_term(A33, B33)\n", "right = right_term(A33, B33)\n", "print(f\"L2 error (3,3)-doubleforms via formula: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL\n" ] }, { "cell_type": "markdown", "id": "ba93ba1f", "metadata": {}, "source": [ "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\n", "\n", "$$\n", "\\begin{align*}\n", "&\\mathcal{S}\\sigma = \\sigma^T-g\\,\\mathrm{Tr}\\sigma,\\\\\n", "&\\mathcal{J}\\sigma=\\sigma^T-\\frac{1}{2}g\\,\\mathrm{Tr}\\sigma.\n", "\\end{align*}\n", "$$\n", "\n", "For all $\\sigma\\in\\Lambda^{1,1}(\\Omega)$ there holds\n", "\n", "$$\n", "\\begin{align*}\n", "&\\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),\\\\\n", "&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).\n", "\\end{align*}\n", "$$\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "50ee3f25", "metadata": {}, "outputs": [], "source": [ "N = 3\n", "\n", "left = mf.S(B11)\n", "right = -1 / math.factorial(N - 2) * mf.inv_star(dg.Wedge(mf.G, B11))\n", "print(f\"L2 error B operator (1,1)-doubleforms: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL\n", "\n", "left = dg.Wedge(mf.G, mf.J(B11))\n", "right = -1 / math.factorial(N - 3) * mf.inv_star(B11)\n", "print(f\"L2 error J operator (1,1)-doubleforms: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL" ] }, { "cell_type": "markdown", "id": "32717c98", "metadata": {}, "source": [ "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\n", "\n", "$$\n", "\\begin{align*}\n", "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).\n", "\\end{align*}\n", "$$\n", "\n", "The maps $\\mathrm{Tr}$ and $s$ commute.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "751e359b", "metadata": {}, "outputs": [], "source": [ "left = mf.Trace(mf.s(A22))\n", "right = mf.s(mf.Trace(A22))\n", "print(\n", " f\"L2 error trace-s operator (2,2)-doubleforms: {l2_error(left, right, mesh3):.6f}\"\n", ")\n", "assert l2_error(left, right, mesh3) < TOL\n", "\n", "left = mf.Trace(mf.s(A33))\n", "right = mf.s(mf.Trace(A33))\n", "print(\n", " f\"L2 error trace-s operator (3,3)-doubleforms: {l2_error(left, right, mesh3):.6f}\"\n", ")\n", "assert l2_error(left, right, mesh3) < TOL" ] }, { "cell_type": "markdown", "id": "3a3fa0f8", "metadata": {}, "source": [ "There holds the Leibniz rule for $s$ for all $\\varphi\\in\\Lambda^{p,q}(\\Omega)$ and $\\Psi\\in\\Lambda^{l,k}(\\Omega)$\n", "\n", "$$\n", "\\begin{align*}\n", "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).\n", "\\end{align*}\n", "$$\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "2a66aab5", "metadata": {}, "outputs": [], "source": [ "def left_term(phi, psi):\n", " return mf.s(dg.Wedge(phi, psi))\n", "\n", "\n", "def right_term(phi, psi):\n", " p = phi.degree_left\n", " q = phi.degree_right\n", " return dg.Wedge(mf.s(phi), psi) + (-1) ** (p + q) * dg.Wedge(phi, mf.s(psi))\n", "\n", "\n", "left = left_term(B11, C11)\n", "right = right_term(B11, C11)\n", "print(f\"L2 error (1,1)-(1,1)-doubleforms: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL\n", "\n", "left = left_term(A22, B12)\n", "right = right_term(A22, B12)\n", "print(f\"L2 error (2,2)-(1,2)-doubleforms: {l2_error(left, right, mesh3):.6f}\")\n", "assert l2_error(left, right, mesh3) < TOL" ] }, { "cell_type": "code", "execution_count": null, "id": "79085a8e", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "ngsolve-dev (3.12.3)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.3" } }, "nbformat": 4, "nbformat_minor": 5 }