{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "69347fef",
   "metadata": {},
   "source": [
    "# Alternating k-forms, wedge product, Hodge star, and exterior (co-) derivative\n",
    "\n",
    "In this notebook, we introduce important objects and identities of exterior calculus.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df2e499a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve import *\n",
    "from netgen.occ import unit_square, unit_cube\n",
    "import ngsdiffgeo as dg\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": "markdown",
   "id": "351b0adb",
   "metadata": {},
   "source": [
    "## 3D Euclidean case\n",
    "We focus on the more interesting case of $N=3$. First we consider the Euclidean manifold, $\\Omega\\subset \\mathbb{R}^3$, with the usual Euclidean inner products. Although the wedge product and exterior derivative are independent of the metric, for the Hodge star and exterior co-derivative this is essential. Later, we consider a non-Euclidean Riemannian manifold with a metric tensor $g\\neq I$.\n",
    "\n",
    "We define the 1-forms $dx,dy,dz\\in \\Lambda^1(\\Omega)$, 2-forms $dx\\wedge dy, dx\\wedge dz, dy\\wedge dz\\in \\Lambda^2(\\Omega)$, and 3-form $dx\\wedge dy\\wedge dz\\in \\Lambda^3(\\Omega)$. The wedge product anti-symmetrizes the dyadic product of the forms. For example there holds\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\alpha\\wedge \\beta = \\alpha \\otimes \\beta - \\beta\\otimes \\alpha,\\qquad \\alpha,\\beta\\in\\Lambda^1(\\Omega).\n",
    "\\end{align*}\n",
    "$$\n",
    "For the wedge product, there holds the associative rule $(\\alpha\\wedge\\beta)\\wedge \\gamma = \\alpha\\wedge(\\beta\\wedge\\gamma)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c044b11a",
   "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",
    "dx_dy = dg.Wedge(dx, dy)\n",
    "dx_dz = dg.Wedge(dx, dz)\n",
    "dy_dz = dg.Wedge(dy, dz)\n",
    "\n",
    "dx_dy2 = dg.TwoForm(CF((0, 1, 0, -1, 0, 0, 0, 0, 0), dims=(3, 3)), dim=3)\n",
    "dx_dz2 = dg.TwoForm(CF((0, 0, 1, 0, 0, 0, -1, 0, 0), dims=(3, 3)), dim=3)\n",
    "dy_dz2 = dg.TwoForm(CF((0, 0, 0, 0, 0, 1, 0, -1, 0), dims=(3, 3)), dim=3)\n",
    "\n",
    "print(f\"L2 error dx_dy - dx_dy2 = {l2_error(dx_dy, dx_dy2, mesh3):.6f}\")\n",
    "assert l2_error(dx_dy, dx_dy2, mesh3) < TOL\n",
    "print(f\"L2 error dx_dz - dx_dz2 = {l2_error(dx_dz, dx_dz2, mesh3):.6f}\")\n",
    "assert l2_error(dx_dz, dx_dz2, mesh3) < TOL\n",
    "print(f\"L2 error dy_dz - dy_dz2 = {l2_error(dy_dz, dy_dz2, mesh3):.6f}\")\n",
    "assert l2_error(dy_dz, dy_dz2, mesh3) < TOL\n",
    "\n",
    "dx_dy_dz = dg.Wedge(dx_dy, dz)\n",
    "dx_dy_dz_as = dg.Wedge(dx, dy_dz)\n",
    "dx_dy_dz2 = dg.ThreeForm(\n",
    "    CF(\n",
    "        (\n",
    "            0,\n",
    "            0,\n",
    "            0,\n",
    "            0,\n",
    "            0,\n",
    "            1,\n",
    "            0,\n",
    "            -1,\n",
    "            0,\n",
    "            0,\n",
    "            0,\n",
    "            -1,\n",
    "            0,\n",
    "            0,\n",
    "            0,\n",
    "            1,\n",
    "            0,\n",
    "            0,\n",
    "            0,\n",
    "            1,\n",
    "            0,\n",
    "            -1,\n",
    "            0,\n",
    "            0,\n",
    "            0,\n",
    "            0,\n",
    "            0,\n",
    "        ),\n",
    "        dims=(3, 3, 3),\n",
    "    ),\n",
    "    dim=3,\n",
    ")\n",
    "print(f\"L2 error dx_dy_dz - dx_dy_dz2 = {l2_error(dx_dy_dz, dx_dy_dz2, mesh3):.6f}\")\n",
    "assert l2_error(dx_dy_dz, dx_dy_dz2, mesh3) < TOL\n",
    "print(f\"L2 error dx_dy_dz - dx_dy_dz_as = {l2_error(dx_dy_dz, dx_dy_dz_as, mesh3):.6f}\")\n",
    "assert l2_error(dx_dy_dz, dx_dy_dz_as, mesh3) < TOL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc3e1a0f",
   "metadata": {},
   "source": [
    "Further, the wedge product is anti-commutative depending on the degree of the involved forms\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\alpha\\wedge\\beta = (-1)^{lk}\\beta\\wedge\\alpha,\\qquad \\alpha\\in\\Lambda^l(\\Omega), \\beta\\in\\Lambda^k(\\Omega),\n",
    "\\end{align*}\n",
    "$$\n",
    "such that $\\alpha\\wedge\\alpha=0$ for all $\\alpha\\in\\Lambda^k(\\Omega)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1275721",
   "metadata": {},
   "outputs": [],
   "source": [
    "dy_dx = dg.Wedge(dy, dx)\n",
    "print(f\"L2 error dy^dx + dx^dy = {l2_norm(dy_dx + dx_dy, mesh3):.6f}\")\n",
    "assert l2_norm(dy_dx + dx_dy, mesh3) < TOL\n",
    "\n",
    "dz_dx_dy = dg.Wedge(dz, dx_dy)\n",
    "print(f\"L2 error dz^(dx^dy) - dx^dy^dz = {l2_norm(dz_dx_dy - dx_dy_dz, mesh3):.6f}\")\n",
    "assert l2_norm(dz_dx_dy - dx_dy_dz, mesh3) < TOL\n",
    "\n",
    "alpha = dg.OneForm(CF((0.3 * x * y, z**2, -0.1 * x)))\n",
    "alpha_wedge_alpha = dg.Wedge(alpha, alpha)\n",
    "print(f\"L2 norm alpha^alpha = {l2_norm(alpha_wedge_alpha, mesh3):.6f}\")\n",
    "assert l2_norm(alpha_wedge_alpha, mesh3) < TOL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f870a4c",
   "metadata": {},
   "source": [
    "Let $X\\in\\mathfrak{X}(\\Omega)$ be a vector field. The contraction (also denoted as interior product)  $\\iota_X:\\Lambda^k(\\Omega)\\to\\Lambda^{k-1}(\\Omega)$, sometimes denoted by $X\\lrcorner$, is defined by\n",
    "$$\n",
    "\\begin{align*}\n",
    "(\\iota_X\\alpha)(Y_1,\\dots,Y_{k-1}) = \\alpha(X,Y_1,\\dots, Y_{k-1}), \\qquad\\alpha\\in\\Lambda^k(\\Omega),\\, Y_i\\in \\mathfrak{X}(\\Omega).\n",
    "\\end{align*}\n",
    "$$ \n",
    "There holds the following identity for the wedge product\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\iota_X(\\beta\\wedge\\gamma) = (\\iota_X\\beta)\\wedge \\gamma + (-1)^k\\beta\\wedge(\\iota_X\\gamma),\\qquad\\beta\\in\\Lambda^k(\\Omega),\\gamma\\in\\Lambda^l(\\Omega).\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c4986bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "mf = dg.RiemannianManifold(Id(3))\n",
    "\n",
    "X = dg.VectorField(CF((0.3 * y * z - x * y, x**2 * z + 0.34 * y**3, -x * y * z)))\n",
    "alpha = dg.OneForm(CF((0.3 * x * y, z**2, -0.1 * x)))\n",
    "beta = 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",
    "iota_X_beta_gamma = mf.Contraction(X, dg.Wedge(alpha, beta))\n",
    "iota_beta_iota_gamma = dg.Wedge(mf.Contraction(X, alpha), beta) - dg.Wedge(\n",
    "    alpha, mf.Contraction(X, beta)\n",
    ")\n",
    "print(\n",
    "    f\"L2 error i_X(alpha^beta) - i_X(alpha)^beta + alpha^i_X(beta) = {l2_error(iota_X_beta_gamma, iota_beta_iota_gamma, mesh3):.6f}\"\n",
    ")\n",
    "assert l2_error(iota_X_beta_gamma, iota_beta_iota_gamma, mesh3) < TOL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "990f32c3",
   "metadata": {},
   "source": [
    "The Hodge star operator $\\star:\\Lambda^{k}(\\Omega)\\to \\Lambda^{N-k}(\\Omega)$ is defined by the identity\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\alpha\\wedge (\\star\\beta) = \\frac{1}{k!}\\langle \\alpha, \\beta\\rangle\\,\\omega,\\qquad \\alpha,\\beta\\in \\Lambda^{k}(\\Omega),\n",
    "\\end{align*}\n",
    "$$\n",
    "where $\\omega$ is the volume form ($dx\\wedge dy\\wedge dz$ in the Euclidean case) and $\\langle \\alpha, \\beta\\rangle$ the inner product of $k$-tensors. Note that the normalization factor $1/k!$ is necessary as in the inner product we sum over all indices. In the literature the inner product of $k$-forms, sometimes denoted by $g^{1}(\\alpha,\\beta)$ is defined by\n",
    "$$\n",
    "\\begin{align*}\n",
    "g^{-1}(\\alpha, \\beta)=\\frac{1}{k!} g^{i_1 j_1}\\cdots g^{i_k j_k}\\, \\alpha_{i_1\\dots i_k}\\, \\beta_{j_1\\dots j_k},\\qquad \\alpha,\\beta\\in \\Lambda^{k}(\\Omega),\n",
    "\\end{align*}\n",
    "$$\n",
    "without the normalization factor, so there holds $\\langle\\alpha,\\beta\\rangle = k!\\,g^{-1}(\\alpha,\\beta)$ for $\\alpha,\\beta\\in\\Lambda^k(\\Omega)$. Applying the Hodge star twice yields the identity up to a sign\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\star\\star \\alpha = (-1)^{k(N-k)}\\alpha, \\alpha\\in\\Lambda^k(\\Omega).\n",
    "\\end{align*}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9087ab48",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Euclidean manifold, metric tensor g=Id\n",
    "mf = dg.RiemannianManifold(Id(3))\n",
    "\n",
    "one = dg.ScalarField(CF(1), dim=3)\n",
    "\n",
    "star_one = mf.star(one)\n",
    "star_dx = mf.star(dx)\n",
    "star_dy = mf.star(dy)\n",
    "star_dz = mf.star(dz)\n",
    "star_dx_dy = mf.star(dx_dy)\n",
    "star_dx_dz = mf.star(dx_dz)\n",
    "star_dy_dz = mf.star(dy_dz)\n",
    "star_dx_dy_dz = mf.star(dx_dy_dz)\n",
    "\n",
    "print(f\"L2 error star(1) - dx^dy^dz = {l2_error(star_one, dx_dy_dz, mesh3):.6f}\")\n",
    "assert l2_error(star_one, dx_dy_dz, mesh3) < TOL\n",
    "print(f\"L2 error star(dx) - dy^dz = {l2_error(star_dx, dy_dz, mesh3):.6f}\")\n",
    "assert l2_error(star_dx, dy_dz, mesh3) < TOL\n",
    "print(f\"L2 error star(dy) - -dz^dx = {l2_error(star_dy, -dx_dz, mesh3):.6f}\")\n",
    "assert l2_error(star_dy, -dx_dz, mesh3) < TOL\n",
    "print(f\"L2 error star(dz) - dx^dy = {l2_error(star_dz, dx_dy, mesh3):.6f}\")\n",
    "assert l2_error(star_dz, dx_dy, mesh3) < TOL\n",
    "print(f\"L2 error star(dx^dy) - dz = {l2_error(star_dx_dy, dz, mesh3):.6f}\")\n",
    "assert l2_error(star_dx_dy, dz, mesh3) < TOL\n",
    "print(f\"L2 error star(dx^dz) - -dy = {l2_error(star_dx_dz, -dy, mesh3):.6f}\")\n",
    "assert l2_error(star_dx_dz, -dy, mesh3) < TOL\n",
    "print(f\"L2 error star(dy^dz) - dx = {l2_error(star_dy_dz, dx, mesh3):.6f}\")\n",
    "assert l2_error(star_dy_dz, dx, mesh3) < TOL\n",
    "print(f\"L2 error star(dx^dy^dz) - 1 = {l2_error(star_dx_dy_dz, one, mesh3):.6f}\")\n",
    "assert l2_error(star_dx_dy_dz, one, mesh3) < TOL\n",
    "\n",
    "alpha = dg.OneForm(CF((0.3 * x * y, z**2, -0.1 * x)))\n",
    "beta = 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",
    "gamma = (0.1 + x * y - 0.2 * x * y * z) * dx_dy_dz\n",
    "\n",
    "star_star_alpha = mf.star(mf.star(alpha))\n",
    "expected_star_star_alpha = alpha\n",
    "star_star_beta = mf.star(mf.star(beta))\n",
    "expected_star_star_beta = beta\n",
    "star_star_gamma = mf.star(mf.star(gamma))\n",
    "expected_star_star_gamma = gamma\n",
    "\n",
    "print(\n",
    "    f\"L2 error star(star(alpha)) - alpha = {l2_error(star_star_alpha, expected_star_star_alpha, mesh3):.6f}\"\n",
    ")\n",
    "assert l2_error(star_star_alpha, expected_star_star_alpha, mesh3) < TOL\n",
    "print(\n",
    "    f\"L2 error star(star(beta)) - beta = {l2_error(star_star_beta, expected_star_star_beta, mesh3):.6f}\",\n",
    ")\n",
    "assert l2_error(star_star_beta, expected_star_star_beta, mesh3) < TOL\n",
    "print(\n",
    "    f\"L2 error star(star(gamma)) - gamma = {l2_error(star_star_gamma, expected_star_star_gamma, mesh3):.6f}\",\n",
    ")\n",
    "assert l2_error(star_star_gamma, expected_star_star_gamma, mesh3) < TOL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b737944",
   "metadata": {},
   "source": [
    "The exterior derivative $d:\\Lambda^k(\\Omega)\\to\\Lambda^{k+1}(\\Omega)$ is defined by its action on vector fields $X_0,\\dots,X_k\\in\\mathfrak{X}(\\Omega)$\n",
    "$$\n",
    "\\begin{align*}\n",
    "(d\\alpha)(X_0,\\dots,X_k) = \\sum_{i=0}^k (-1)^i X_i\\big(\\alpha(X_0,\\dots,\\hat{X_i},\\dots,X_k)\\big)\n",
    " + \\sum_{0\\leq i<j\\leq k} (-1)^{i+j} \\alpha([X_i,X_j],X_0,\\dots,\\hat{X_i},\\dots,\\hat{X_j},\\dots,X_k),\n",
    "\\end{align*}\n",
    "$$\n",
    "which in coordinates reduces to the alternating sum of partial derivatives of the coefficients. Here, $[X,Y]$ denotes the Lie-bracket (or commutator), which is defined by the identity $[X,Y]f = X(Y(f))-Y(X(f))$ for all functions $f$. Further, $\\hat\n",
    "X$ means that this argument gets skipped.\n",
    "\n",
    " Applying the exterior derivative twice yields zero, $dd \\alpha = 0$ for $\\alpha\\in\\Lambda^k(\\Omega)$, i.e., $d$ is nilpotent. Further, there holds the following Leibnitz rule with the wedge product\n",
    "$$\n",
    "\\begin{align*}\n",
    "d(\\alpha\\wedge\\beta) = d\\alpha\\wedge\\beta + (-1)^k\\alpha\\wedge d\\beta,\\quad\\alpha\\in\\Lambda^k(\\Omega),\\beta\\in\\Lambda^l(\\Omega).\n",
    "\\end{align*}\n",
    "$$\n",
    "For 1-forms there also holds $d \\alpha = \\nabla \\alpha - (\\nabla\\alpha)^T$ with $\\nabla$ denoting the Levi-Civita connection (covariant derivative) as the involved Christoffel symbols cancel out due to their symmetry.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8f4c8ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "g = dg.ScalarField(x * y * z, dim=3)\n",
    "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",
    "\n",
    "dg_g = dg.d(g)\n",
    "expected_dg_g = dg.OneForm(CF((y * z, x * z, x * y)))\n",
    "print(f\"|d g - expected||_L2: {l2_norm(dg_g - expected_dg_g, mesh3):.6f}\")\n",
    "assert l2_norm(dg_g - expected_dg_g, mesh3) < TOL\n",
    "ddg_g = dg.d(dg_g)\n",
    "print(f\"d^2 g = {l2_norm(ddg_g, mesh3):.6f}\")\n",
    "assert l2_norm(ddg_g, mesh3) < TOL\n",
    "\n",
    "dalpha = dg.d(alpha)\n",
    "expected_dalpha = mf.CovDeriv(alpha) - mf.CovDeriv(alpha).trans\n",
    "print(f\"||d alpha - expected||_L2: {l2_norm(dalpha - expected_dalpha, mesh3):.6f}\")\n",
    "assert l2_norm(dalpha - expected_dalpha, mesh3) < TOL\n",
    "ddalpha = dg.d(dalpha)\n",
    "print(f\"d^2 alpha = {l2_norm(ddalpha, mesh3):.6f}\")\n",
    "assert l2_norm(ddalpha, mesh3) < TOL\n",
    "\n",
    "d_a_w_b = dg.d(dg.Wedge(alpha, beta))\n",
    "expected_d_a_w_b = dg.Wedge(dalpha, beta) - dg.Wedge(alpha, dg.d(beta))\n",
    "print(\n",
    "    f\"||d(alpha^beta) - expected||_L2: {l2_norm(d_a_w_b - expected_d_a_w_b, mesh3):.6f}\"\n",
    ")\n",
    "assert l2_norm(d_a_w_b - expected_d_a_w_b, mesh3) < TOL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98ea0eac",
   "metadata": {},
   "source": [
    "For a vector field $X\\in\\mathfrak{X}(\\Omega)$ the covariant divergence and curl (in 3D) can be written in terms of the exterior derivative, Hodge star, and the musical isomorphisms relating 1-forms with vector fields\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\mathrm{div}\\, X = \\star d\\star (X^\\flat),\\qquad \\mathrm{curl}\\,X = (\\star d(X^\\flat))^\\sharp.\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "460f57b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "mf = dg.RiemannianManifold(Id(3))\n",
    "X = dg.VectorField(CF((0.3 * y * z - x * y, x**2 * z + 0.34 * y**3, -x * y * z)))\n",
    "\n",
    "div_X = mf.star(dg.d(mf.star(mf.Lower(X))))\n",
    "div_X_expected = mf.CovDiv(X)\n",
    "print(f\"||div X - expected||_L2: {l2_norm(div_X - div_X_expected, mesh3):.6f}\")\n",
    "assert l2_norm(div_X - div_X_expected, mesh3) < TOL\n",
    "\n",
    "curl_X = mf.Raise(mf.star(dg.d(mf.Lower(X))))\n",
    "curl_X_expected = mf.CovCurl(X)\n",
    "print(f\"||curl X - expected||_L2: {l2_norm(curl_X - curl_X_expected, mesh3):.6f}\")\n",
    "assert l2_norm(curl_X - curl_X_expected, mesh3) < TOL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e48b5a3",
   "metadata": {},
   "source": [
    "Further, there holds the Stokes theorem for $\\alpha\\in \\Lambda^{N-1}(\\Omega)$\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\int_{\\Omega}d\\alpha = \\int_{\\partial\\Omega}\\iota^*\\alpha,\n",
    "\\end{align*}\n",
    "$$\n",
    "where $\\iota^*$ is the pullback of the inclusion $\\iota:\\partial\\Omega\\to\\Omega$. We use that $d\\alpha = \\star(d\\alpha)\\omega$  and $\\iota^*\\alpha = \\star_\\partial(\\iota^*\\alpha)\\omega_{\\partial}$ for $\\alpha\\in \\Lambda^{N-1}(\\Omega)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c39f88d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import ngsolve\n",
    "\n",
    "mf = dg.RiemannianManifold(Id(3))\n",
    "beta = 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",
    "            0.2 * x * z,\n",
    "            0,\n",
    "            -0.2 * x * z,\n",
    "            0,\n",
    "        ),\n",
    "        dims=(3, 3),\n",
    "    )\n",
    ")\n",
    "\n",
    "\n",
    "left = Integrate(mf.star(dg.d(beta)) * ngsolve.dx(bonus_intorder=5), mesh3)\n",
    "right = Integrate(\n",
    "    mf.star(beta, BND) * ngsolve.dx(bonus_intorder=5, element_boundary=True), mesh3\n",
    ")\n",
    "\n",
    "print(f\"Left integral: {left:.6f}, right integral: {right:.6f}\")\n",
    "assert abs(left - right) < TOL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc96a941",
   "metadata": {},
   "source": [
    "The exterior coderivative $\\delta:\\Lambda^k(\\Omega)\\to\\Lambda^{k-1}(\\Omega)$ is defined as the adjoint of the exterior derivative with respect to the $L^2$ inner product, i.e., for all $\\alpha\\in\\Lambda^k(\\Omega)$ and $\\beta\\in \\Lambda^{k+1}(\\Omega)$ with compact support there holds\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\int_{\\Omega}\\frac{1}{(k+1)!}\\langle d\\alpha,\\beta\\rangle\\,\\omega = \\int_{\\Omega}\\frac{1}{k!}\\langle \\alpha,\\delta\\beta\\rangle\\,\\omega,\n",
    "\\end{align*}\n",
    "$$\n",
    "where $\\omega$ and $\\langle\\cdot,\\cdot\\rangle$ are again the volume form and the inner product for $k$-tensor fields. There holds the identity $\\delta = (-1)^{N(k+1)+1}\\,\\star d \\star$. From it, there follows that $\\delta$ is also nilpotent, $\\delta\\delta=0$.\n",
    "\n",
    " Using Stokes' theorem together with  yields the following integration-by-parts formula that holds for all $\\alpha\\in\\Lambda^k(\\Omega)$ and $\\beta\\in \\Lambda^{k+1}(\\Omega)$\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\int_{\\Omega}\\frac{1}{(k+1)!}\\langle d\\alpha,\\beta\\rangle\\,\\omega = \\int_{\\Omega}\\frac{1}{k!}\\langle \\alpha,\\delta\\beta\\rangle\\,\\omega + \\int_{\\partial \\Omega}\\frac{1}{k!}\\langle \\alpha,\\iota_n\\beta\\rangle\\,\\omega_{\\partial},\n",
    "\\end{align*}\n",
    "$$\n",
    "where $\\omega_{\\partial}$ is the volume form on the boundary, $n$ the outer unit normal vector, and for a vector field $X\\in\\mathfrak{X}(\\Omega)$ the contraction is $\\iota_X$. Further, the inner product in the boundary term has to be understood with respect to the tangential metric $g_F$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64a20e9d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import ngsolve\n",
    "\n",
    "mf = dg.RiemannianManifold(Id(3))\n",
    "\n",
    "f = dg.ScalarField(200 * x * (1 - x) * y * (1 - y) * z * (1 - z), dim=3)\n",
    "alpha = dg.OneForm(CF((0.3 * x * y, z**2, -0.1 * x)))\n",
    "\n",
    "left = Integrate(\n",
    "    mf.InnerProduct(alpha, dg.d(f)) * mf.VolumeForm(VOL) * ngsolve.dx(bonus_intorder=5),\n",
    "    mesh3,\n",
    ")\n",
    "right = Integrate(\n",
    "    mf.InnerProduct(f, mf.delta(alpha))\n",
    "    * mf.VolumeForm(VOL)\n",
    "    * ngsolve.dx(bonus_intorder=5),\n",
    "    mesh3,\n",
    ")\n",
    "\n",
    "print(f\"Compact support: left: {left:.6f}, right: {right:.6f}\")\n",
    "assert abs(left - right) < TOL\n",
    "\n",
    "f = dg.ScalarField(20 * x * y * (1 - y) * (1 - z), dim=3)\n",
    "alpha = dg.OneForm(CF((0.3 * x * y, z**2, -0.1 * x)))\n",
    "\n",
    "left = Integrate(\n",
    "    mf.InnerProduct(alpha, dg.d(f)) * mf.VolumeForm(VOL) * ngsolve.dx(bonus_intorder=5),\n",
    "    mesh3,\n",
    ")\n",
    "right = Integrate(\n",
    "    mf.InnerProduct(f, mf.delta(alpha))\n",
    "    * mf.VolumeForm(VOL)\n",
    "    * ngsolve.dx(bonus_intorder=5),\n",
    "    mesh3,\n",
    ") + Integrate(\n",
    "    f\n",
    "    * mf.Contraction(mf.normal, alpha)\n",
    "    * mf.VolumeForm(BND)\n",
    "    * ngsolve.dx(element_boundary=True, bonus_intorder=5),\n",
    "    mesh3,\n",
    ")\n",
    "print(f\"Integration by parts: left: {left:.6f}, right: {right:.6f}\")\n",
    "assert abs(left - right) < TOL\n",
    "\n",
    "alpha = dg.OneForm(CF((0.3 * x * y, z**2, -0.1 * x)))\n",
    "beta = 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",
    "left = Integrate(\n",
    "    0.5\n",
    "    * mf.InnerProduct(beta, dg.d(alpha))\n",
    "    * mf.VolumeForm(VOL)\n",
    "    * ngsolve.dx(bonus_intorder=5),\n",
    "    mesh3,\n",
    ")\n",
    "right = Integrate(\n",
    "    mf.InnerProduct(alpha, mf.delta(beta))\n",
    "    * mf.VolumeForm(VOL)\n",
    "    * ngsolve.dx(bonus_intorder=5),\n",
    "    mesh3,\n",
    ") + Integrate(\n",
    "    mf.InnerProduct(alpha, mf.Contraction(mf.normal, beta), BND)\n",
    "    * mf.VolumeForm(BND)\n",
    "    * ngsolve.dx(element_boundary=True, bonus_intorder=5),\n",
    "    mesh3,\n",
    ")\n",
    "print(f\"Integration by parts: left: {left:.6f}, right: {right:.6f}\")\n",
    "assert abs(left - right) < TOL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a36c6469",
   "metadata": {},
   "source": [
    "The Laplace-Beltrami operator can be expressed by the exterior (co-)derivative \n",
    "$$\n",
    "\\begin{align*}\n",
    "\\Delta \\alpha = -(d\\delta +\\delta d)\\alpha,\\qquad \\alpha\\in \\Lambda^k(\\Omega).\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2623c58",
   "metadata": {},
   "outputs": [],
   "source": [
    "mf = dg.RiemannianManifold(Id(3))\n",
    "\n",
    "f = dg.ScalarField(20 * x * y * (1 - y) * (1 - z), dim=3)\n",
    "alpha = dg.OneForm(CF((0.3 * x * y, z**2, -0.1 * x)))\n",
    "\n",
    "laplace_beltrami = -(dg.d(mf.delta(f)) + mf.delta(dg.d(f)))\n",
    "expected = mf.CovLaplace(f)\n",
    "print(f\"Delta f: {l2_error(laplace_beltrami, expected, mesh3):.6f}\")\n",
    "assert l2_error(laplace_beltrami, expected, mesh3) < TOL\n",
    "\n",
    "laplace_beltrami = -(dg.d(mf.delta(alpha)) + mf.delta(dg.d(alpha)))\n",
    "expected = mf.CovLaplace(alpha)\n",
    "print(f\"Delta alpha: {l2_error(laplace_beltrami, expected, mesh3):.6f}\")\n",
    "assert l2_error(laplace_beltrami, expected, mesh3) < TOL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e7e6c9f",
   "metadata": {},
   "source": [
    "## 3D with general Riemannian manifold\n",
    "Nest, we consider a general Riemannian manifold. Then the Hodge star and exterior co-derivative read in coordinates\n",
    "$$\n",
    "\\begin{align*}\n",
    ".\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "Adapting the volume forms, $\\omega=\\sqrt{\\det g}\\,dx\\wedge dy\\wedge dz$, $\\omega_{\\partial}=\\sqrt{\\|\\mathrm{cof}(g)\\|_F}\\,dt_1\\wedge dt_2$, the $g$ inner product $\\langle\\cdot,\\cdot\\rangle= g(\\cdot,\\cdot)$ extended from vector fields to general tensor fields, and the $g$-normal vector $n=n_g$ the above identities hold analogously."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4d1fe16",
   "metadata": {},
   "outputs": [],
   "source": [
    "import ngsolve\n",
    "\n",
    "\n",
    "mf = dg.RiemannianManifold(dg.WarpedProduct().metric)\n",
    "\n",
    "################ covariant div and curl ################\n",
    "X = dg.VectorField(CF((0.3 * y * z - x * y, x**2 * z + 0.34 * y**3, -x * y * z)))\n",
    "\n",
    "div_X = mf.star(dg.d(mf.star(mf.Lower(X))))\n",
    "div_X_expected = mf.CovDiv(X)\n",
    "print(f\"||div X - expected||_L2: {l2_norm(div_X - div_X_expected, mesh3):.6f}\")\n",
    "assert l2_norm(div_X - div_X_expected, mesh3) < TOL\n",
    "\n",
    "curl_X = mf.Raise(mf.star(dg.d(mf.Lower(X))))\n",
    "curl_X_expected = mf.CovCurl(X)\n",
    "print(f\"||curl X - expected||_L2: {l2_norm(curl_X - curl_X_expected, mesh3):.6f}\")\n",
    "assert l2_norm(curl_X - curl_X_expected, mesh3) < TOL\n",
    "\n",
    "################ exterior co-derivative and integration-by-parts ################\n",
    "f = dg.ScalarField(200 * x * (1 - x) * y * (1 - y) * z * (1 - z), dim=3)\n",
    "alpha = dg.OneForm(CF((0.3 * x * y, z**2, -0.1 * x)))\n",
    "\n",
    "left = Integrate(\n",
    "    mf.InnerProduct(alpha, dg.d(f)) * mf.VolumeForm(VOL) * ngsolve.dx(bonus_intorder=5),\n",
    "    mesh3,\n",
    ")\n",
    "right = Integrate(\n",
    "    mf.InnerProduct(f, mf.delta(alpha))\n",
    "    * mf.VolumeForm(VOL)\n",
    "    * ngsolve.dx(bonus_intorder=5),\n",
    "    mesh3,\n",
    ")\n",
    "print(f\"Compact support: left: {left:.6f}, right: {right:.6f}\")\n",
    "assert abs(left - right) < TOL\n",
    "\n",
    "f = dg.ScalarField(20 * x * y * (1 - y) * (1 - z), dim=3)\n",
    "alpha = dg.OneForm(CF((0.3 * x * y, z**2, -0.1 * x)))\n",
    "\n",
    "left = Integrate(\n",
    "    mf.InnerProduct(alpha, dg.d(f)) * mf.VolumeForm(VOL) * ngsolve.dx(bonus_intorder=5),\n",
    "    mesh3,\n",
    ")\n",
    "right = Integrate(\n",
    "    mf.InnerProduct(f, mf.delta(alpha))\n",
    "    * mf.VolumeForm(VOL)\n",
    "    * ngsolve.dx(bonus_intorder=5),\n",
    "    mesh3,\n",
    ") + Integrate(\n",
    "    f\n",
    "    * mf.Contraction(mf.normal, alpha)\n",
    "    * mf.VolumeForm(BND)\n",
    "    * ngsolve.dx(element_boundary=True, bonus_intorder=5),\n",
    "    mesh3,\n",
    ")\n",
    "print(f\"Integration by parts: left: {left:.6f}, right: {right:.6f}\")\n",
    "assert abs(left - right) < TOL\n",
    "\n",
    "alpha = dg.OneForm(CF((0.3 * x * y, z**2, -0.1 * x)))\n",
    "beta = 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",
    "left = Integrate(\n",
    "    0.5\n",
    "    * mf.InnerProduct(beta, dg.d(alpha))\n",
    "    * mf.VolumeForm(VOL)\n",
    "    * ngsolve.dx(bonus_intorder=5),\n",
    "    mesh3,\n",
    ")\n",
    "right = Integrate(\n",
    "    mf.InnerProduct(alpha, mf.delta(beta))\n",
    "    * mf.VolumeForm(VOL)\n",
    "    * ngsolve.dx(bonus_intorder=5),\n",
    "    mesh3,\n",
    ") + Integrate(\n",
    "    mf.InnerProduct(alpha, mf.Contraction(mf.normal, beta), BND)\n",
    "    * mf.VolumeForm(BND)\n",
    "    * ngsolve.dx(element_boundary=True, bonus_intorder=5),\n",
    "    mesh3,\n",
    ")\n",
    "print(f\"Integration by parts: left: {left:.6f}, right: {right:.6f}\")\n",
    "assert abs(left - right) < TOL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab2f4c8e",
   "metadata": {},
   "source": [
    "## 2D Euclidean case\n",
    "\n",
    "We also present some identities in the two-dimensional Euclidean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e091600d",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh2 = Mesh(unit_square.GenerateMesh(maxh=2))\n",
    "\n",
    "dx = dg.OneForm(CF((1, 0)))\n",
    "dy = dg.OneForm(CF((0, 1)))\n",
    "\n",
    "dx_dy = dg.Wedge(dx, dy)\n",
    "\n",
    "dx_dy2 = dg.TwoForm(CF((0, 1, -1, 0), dims=(2, 2)))\n",
    "\n",
    "print(f\"L2 error dx_dy - dx_dy2 = {l2_error(dx_dy, dx_dy2, mesh2):.6f}\")\n",
    "assert l2_error(dx_dy, dx_dy2, mesh2) < TOL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d372cf6",
   "metadata": {},
   "outputs": [],
   "source": [
    "dy_dx = dg.Wedge(dy, dx)\n",
    "print(f\"L2 error dy^dx + dx^dy = {l2_norm(dy_dx + dx_dy, mesh2):.6f}\")\n",
    "assert l2_norm(dy_dx + dx_dy, mesh2) < TOL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c68baae",
   "metadata": {},
   "outputs": [],
   "source": [
    "mf = dg.RiemannianManifold(Id(2))\n",
    "one = dg.ScalarField(CF(1), dim=2)\n",
    "star_one = mf.star(one)\n",
    "star_dx = mf.star(dx)\n",
    "star_dy = mf.star(dy)\n",
    "star_dx_dy = mf.star(dx_dy)\n",
    "print(f\"L2 error star(1) - dx^dy = {l2_error(star_one, dx_dy, mesh2):.6f}\")\n",
    "assert l2_error(star_one, dx_dy, mesh2) < TOL\n",
    "print(f\"L2 error star(dx) - dy = {l2_error(star_dx, dy, mesh2):.6f}\")\n",
    "assert l2_error(star_dx, dy, mesh2) < TOL\n",
    "print(f\"L2 error star(dy) - -dx = {l2_error(star_dy, -dx, mesh2):.6f}\")\n",
    "assert l2_error(star_dy, -dx, mesh2) < TOL\n",
    "print(f\"L2 error star(dx^dy) - 1 = {l2_error(star_dx_dy, one, mesh2):.6f}\")\n",
    "assert l2_error(star_dx_dy, one, mesh2) < TOL\n",
    "\n",
    "left = dg.Wedge(dx, star_dx)\n",
    "right = mf.InnerProduct(dx, dx) * dx_dy\n",
    "print(f\"L2 error dx^star(dx) - <dx,dx>*dx^dy = {l2_error(left, right, mesh2):.6f}\")\n",
    "assert l2_error(left, right, mesh2) < TOL\n",
    "alpha = dg.OneForm(CF((0.3 * x * y, -0.1 * x)))\n",
    "beta = (y - 0.1 * x**2) * dx_dy\n",
    "\n",
    "star_star_alpha = mf.star(mf.star(alpha))\n",
    "expected_star_star_alpha = -alpha\n",
    "star_star_beta = mf.star(mf.star(beta))\n",
    "expected_star_star_beta = beta\n",
    "\n",
    "print(\n",
    "    f\"L2 error star(star(alpha)) + alpha = {l2_error(star_star_alpha, expected_star_star_alpha, mesh2):.6f}\"\n",
    ")\n",
    "assert l2_error(star_star_alpha, expected_star_star_alpha, mesh2) < TOL\n",
    "print(\n",
    "    f\"L2 error star(star(beta)) - beta = {l2_error(star_star_beta, expected_star_star_beta, mesh2):.6f}\"\n",
    ")\n",
    "assert l2_error(star_star_beta, expected_star_star_beta, mesh2) < TOL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01cd3650",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh2 = Mesh(unit_square.GenerateMesh(maxh=0.25))\n",
    "\n",
    "# scalar 0-form and 1-form\n",
    "f = dg.ScalarField(x * y, dim=2)\n",
    "alpha = dg.OneForm(CF((x**2, y**2)))\n",
    "\n",
    "df = dg.d(f)\n",
    "ddf = dg.d(df)\n",
    "print(f\"||d^2 f||_L2: {l2_norm(ddf, mesh2):.6f}\")\n",
    "assert l2_norm(ddf, mesh2) < TOL\n",
    "\n",
    "# Leibniz rule d(f alpha) = df^alpha + f d alpha\n",
    "left = dg.d(dg.Wedge(f, alpha))\n",
    "right = dg.Wedge(df, alpha) + dg.Wedge(f, dg.d(alpha))\n",
    "print(f\"Leibniz residual L2: {l2_error(left, right, mesh2):.6f}\")\n",
    "assert l2_error(left, right, mesh2) < TOL\n",
    "\n",
    "# Alternation of wedge of 1-forms\n",
    "beta = dg.OneForm(CF((y, x)))\n",
    "ab = dg.Wedge(alpha, beta)\n",
    "ba = dg.Wedge(beta, alpha)\n",
    "print(f\"||alpha^beta + beta^alpha||_L2: {l2_norm(ab + ba, mesh2):.6f}\")\n",
    "assert l2_norm(ab + ba, mesh2) < TOL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ngsolve-dev",
   "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
}