{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f4681d7b",
   "metadata": {},
   "source": [
    "# Exterior covariant derivatives of double forms\n",
    "\n",
    "This notebook continues [Double forms](14_double_forms.ipynb). We introduce the exterior covariant derivative and coderivative in either slot, and verify their basic algebraic and integration-by-parts identities.\n",
    "\n",
    "$$\n",
    "\\def\\owedge{\\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}}}\n",
    "\\def\\Tr{\\mathrm{Tr}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f584b8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import ngsolve\n",
    "from ngsolve import (\n",
    "    BND,\n",
    "    VOL,\n",
    "    CF,\n",
    "    Id,\n",
    "    InnerProduct,\n",
    "    Integrate,\n",
    "    Mesh,\n",
    "    TaskManager,\n",
    "    sqrt,\n",
    "    x,\n",
    "    y,\n",
    "    z,\n",
    ")\n",
    "from netgen.occ import unit_cube\n",
    "import ngsdiffgeo as dg\n",
    "\n",
    "TOL = 1e-7\n",
    "mesh = Mesh(unit_cube.GenerateMesh(maxh=2))\n",
    "mf = dg.RiemannianManifold(dg.Heisenberg().metric, normal_sign=-1)\n",
    "\n",
    "\n",
    "def l2_error(a, b, vb=VOL, bonus_intorder=4):\n",
    "    dX = (\n",
    "        ngsolve.dx(element_boundary=True, bonus_intorder=bonus_intorder)\n",
    "        if vb == BND\n",
    "        else ngsolve.dx(bonus_intorder=bonus_intorder)\n",
    "    )\n",
    "    return sqrt(Integrate(InnerProduct(a - b, a - b) * dX, mesh))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71c35ff8",
   "metadata": {},
   "source": [
    "We use smooth double forms on a three-dimensional non-Euclidean manifold. A double form of bidegree $(p,q)$ is represented by `dg.DoubleForm(..., p=p, q=q, dim=3)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d23bd362",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = dg.OneForm(CF((0.3 * x * y, z**2, -0.1 * x)))\n",
    "beta = dg.OneForm(CF((0.3 * z * y, x * z**2, y**2)))\n",
    "gamma = dg.OneForm(CF((y * z - x, x**2 + z, x * y)))\n",
    "eta = dg.OneForm(CF((z + x, y - z, x * y)))\n",
    "\n",
    "beta_gamma = dg.Wedge(beta, gamma)\n",
    "alpha_beta = dg.Wedge(alpha, beta)\n",
    "\n",
    "A11 = dg.DoubleForm(dg.Einsum(\"i,j->ij\", alpha, beta), p=1, q=1, dim=3)\n",
    "B11 = dg.DoubleForm(dg.Einsum(\"i,j->ij\", gamma, eta), p=1, q=1, dim=3)\n",
    "A12 = dg.DoubleForm(dg.Einsum(\"i,jk->ijk\", alpha, beta_gamma), p=1, q=2, dim=3)\n",
    "A21 = dg.DoubleForm(dg.Einsum(\"ij,k->ijk\", alpha_beta, gamma), p=2, q=1, dim=3)\n",
    "A22 = dg.Wedge(A11, B11)\n",
    "X = dg.VectorField(CF((x - 0.1 * y * z, y + 0.2 * z * x, z - 0.3 * x * y)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04d321ae",
   "metadata": {},
   "source": [
    "## Left and right slots\n",
    "\n",
    "For $\\varphi\\in\\Lambda^{p,q}(\\Omega)$, the operators `mf.d_cov` and `mf.delta_cov` act on the left slot by default. Set `slot=\"right\"` to act on the right slot. Thus,\n",
    "\n",
    "$$\n",
    "d^\\nabla:\\Lambda^{p,q}\\to\\Lambda^{p+1,q},\\qquad\n",
    "d^{\\nabla,\\prime}:\\Lambda^{p,q}\\to\\Lambda^{p,q+1},\n",
    "$$\n",
    "\n",
    "and the two coderivatives lower the corresponding degree. Transposition exchanges the slots:\n",
    "\n",
    "$$\n",
    "d^{\\nabla,\\prime}\\varphi=\\bigl(d^\\nabla(\\varphi^T)\\bigr)^T,\\qquad\n",
    "\\delta^{\\nabla,\\prime}\\varphi=\\bigl(\\delta^\\nabla(\\varphi^T)\\bigr)^T.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f7b5830",
   "metadata": {},
   "outputs": [],
   "source": [
    "with TaskManager():\n",
    "    for phi in (A11, A12, A21, A22):\n",
    "        d_right = mf.d_cov(phi, slot=\"right\")\n",
    "        d_right_from_transpose = mf.d_cov(phi.trans, slot=\"left\").trans\n",
    "        delta_right = mf.delta_cov(phi, slot=\"right\")\n",
    "        delta_right_from_transpose = mf.delta_cov(phi.trans, slot=\"left\").trans\n",
    "\n",
    "        assert l2_error(d_right, d_right_from_transpose) < TOL\n",
    "        assert l2_error(delta_right, delta_right_from_transpose) < TOL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aedc1624",
   "metadata": {},
   "source": [
    "## Covariant divergence and graph compilation\n",
    "\n",
    "`mf.CovDiv` contracts the covariant derivative in the selected slot. On a double form, the corresponding covariant coderivative is its negative:\n",
    "\n",
    "$$\n",
    "\\delta^\\nabla\\varphi=-\\operatorname{CovDiv}_{L}\\varphi,\\qquad\n",
    "\\delta^{\\nabla,\\prime}\\varphi=-\\operatorname{CovDiv}_{R}\\varphi.\n",
    "$$\n",
    "\n",
    "For a complicated coefficient expression that is evaluated repeatedly, `compile_inner=\"graph\"` can compile the input graph shared by the gradient and connection terms of `d_cov` and `delta_cov`. This is an opt-in optimization for coefficient expressions. Trial and test proxy expressions should use the default mode."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "403739ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "with TaskManager():\n",
    "    for slot in (\"left\", \"right\"):\n",
    "        assert l2_error(mf.CovDiv(A11, slot=slot), -mf.delta_cov(A11, slot=slot)) < TOL\n",
    "        assert (\n",
    "            l2_error(\n",
    "                mf.d_cov(A11, slot=slot, compile_inner=\"graph\"),\n",
    "                mf.d_cov(A11, slot=slot),\n",
    "            )\n",
    "            < TOL\n",
    "        )\n",
    "        assert (\n",
    "            l2_error(\n",
    "                mf.delta_cov(A11, slot=slot, compile_inner=\"graph\"),\n",
    "                mf.delta_cov(A11, slot=slot),\n",
    "            )\n",
    "            < TOL\n",
    "        )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "392c830e",
   "metadata": {},
   "source": [
    "## Leibniz rules\n",
    "\n",
    "The exterior covariant derivative is a graded derivation in the selected slot. For $\\varphi\\in\\Lambda^{p,q}$ and $\\psi\\in\\Lambda^{k,l}$,\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "d^\\nabla(\\varphi\\owedge\\psi)&=d^\\nabla\\varphi\\owedge\\psi+(-1)^p\\varphi\\owedge d^\\nabla\\psi,\\\\\n",
    "d^{\\nabla,\\prime}(\\varphi\\owedge\\psi)&=d^{\\nabla,\\prime}\\varphi\\owedge\\psi+(-1)^q\\varphi\\owedge d^{\\nabla,\\prime}\\psi.\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79dc4730",
   "metadata": {},
   "outputs": [],
   "source": [
    "def d_leibniz(phi, psi, slot):\n",
    "    degree = phi.degree_left if slot == \"left\" else phi.degree_right\n",
    "    return dg.Wedge(mf.d_cov(phi, slot=slot), psi) + (-1) ** degree * dg.Wedge(\n",
    "        phi, mf.d_cov(psi, slot=slot)\n",
    "    )\n",
    "\n",
    "\n",
    "with TaskManager():\n",
    "    for phi, psi in ((A11, B11), (A12, B11), (A21, B11)):\n",
    "        for slot in (\"left\", \"right\"):\n",
    "            left = mf.d_cov(dg.Wedge(phi, psi), slot=slot)\n",
    "            right = d_leibniz(phi, psi, slot)\n",
    "            assert l2_error(left, right) < TOL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "187a3691",
   "metadata": {},
   "source": [
    "The trace relates the two coderivatives to the derivative in the opposite slot:\n",
    "\n",
    "$$\n",
    "\\delta^{\\nabla,\\prime}=-\\Tr d^\\nabla-d^\\nabla\\Tr,\\qquad\n",
    "\\delta^\\nabla=-\\Tr d^{\\nabla,\\prime}-d^{\\nabla,\\prime}\\Tr.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "089ffae9",
   "metadata": {},
   "outputs": [],
   "source": [
    "with TaskManager():\n",
    "    delta_right = mf.delta_cov(A22, slot=\"right\")\n",
    "    trace_formula_right = -mf.Trace(mf.d_cov(A22, slot=\"left\")) - mf.d_cov(\n",
    "        mf.Trace(A22), slot=\"left\"\n",
    "    )\n",
    "    delta_left = mf.delta_cov(A22, slot=\"left\")\n",
    "    trace_formula_left = -mf.Trace(mf.d_cov(A22, slot=\"right\")) - mf.d_cov(\n",
    "        mf.Trace(A22), slot=\"right\"\n",
    "    )\n",
    "\n",
    "    assert l2_error(delta_right, trace_formula_right) < TOL\n",
    "    assert l2_error(delta_left, trace_formula_left) < TOL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7888e51",
   "metadata": {},
   "source": [
    "The interior product is also a graded derivation. The method `ContractSlot` preserves the unselected slot and contracts the selected slot with a vector field."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c0796e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def contraction_leibniz(phi, psi, slot):\n",
    "    degree = phi.degree_left if slot == \"left\" else phi.degree_right\n",
    "    return dg.Wedge(mf.ContractSlot(phi, X, slot=slot), psi) + (\n",
    "        -1\n",
    "    ) ** degree * dg.Wedge(phi, mf.ContractSlot(psi, X, slot=slot))\n",
    "\n",
    "\n",
    "with TaskManager():\n",
    "    for slot in (\"left\", \"right\"):\n",
    "        left = mf.ContractSlot(dg.Wedge(A11, B11), X, slot=slot)\n",
    "        right = contraction_leibniz(A11, B11, slot)\n",
    "        assert l2_error(left, right) < TOL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd9ca96e",
   "metadata": {},
   "source": [
    "## Integration by parts\n",
    "\n",
    "For $\\chi\\in\\Lambda^{p,q}(\\Omega)$ and $\\psi\\in\\Lambda^{p+1,q}(\\Omega)$, the left-slot operators satisfy\n",
    "\n",
    "$$\n",
    "\\int_\\Omega \\langle d^\\nabla\\chi,\\psi\\rangle\\,\\omega\n",
    "=\\int_\\Omega \\langle\\chi,\\delta^\\nabla\\psi\\rangle\\,\\omega\n",
    "-\\int_{\\partial\\Omega}\\langle P_F\\chi,P_n\\psi\\rangle\\,\\omega_{\\partial\\Omega}.\n",
    "$$\n",
    "\n",
    "Here `ProjectDoubleForm(..., left=\"F\")` projects the left slot tangentially and `left=\"n\"` contracts it with the unit normal. We chose `normal_sign=-1`, so the normal is inward and the boundary term has the displayed minus sign. The form inner product includes the factorial normalization in each alternating slot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "030d0b64",
   "metadata": {},
   "outputs": [],
   "source": [
    "chi = A11\n",
    "psi = A21\n",
    "\n",
    "with TaskManager():\n",
    "    volume_left = Integrate(\n",
    "        mf.InnerProduct(mf.d_cov(chi), psi, forms=True)\n",
    "        * mf.VolumeForm(VOL)\n",
    "        * ngsolve.dx(bonus_intorder=6),\n",
    "        mesh,\n",
    "    )\n",
    "    volume_right = Integrate(\n",
    "        mf.InnerProduct(chi, mf.delta_cov(psi), forms=True)\n",
    "        * mf.VolumeForm(VOL)\n",
    "        * ngsolve.dx(bonus_intorder=6),\n",
    "        mesh,\n",
    "    )\n",
    "    boundary = Integrate(\n",
    "        mf.InnerProduct(\n",
    "            mf.ProjectDoubleForm(chi, left=\"F\"),\n",
    "            mf.ProjectDoubleForm(psi, left=\"n\"),\n",
    "            forms=True,\n",
    "        )\n",
    "        * mf.VolumeForm(BND)\n",
    "        * ngsolve.dx(element_boundary=True, bonus_intorder=6),\n",
    "        mesh,\n",
    "    )\n",
    "\n",
    "assert abs(volume_left - volume_right + boundary) < 10 * TOL\n",
    "print(f\"integration-by-parts error: {abs(volume_left - volume_right + boundary):.3e}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 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.13.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}