{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c166542e",
   "metadata": {},
   "source": [
    "# Differential operators and curvature on hypersurfaces\n",
    "\n",
    "This notebook studies the geometry induced on a boundary facet $F$ of a three-dimensional Riemannian manifold $(\\Omega,g)$. We compare ambient and tangential metrics and connections, verify the Gauss formula and Codazzi equation, and use boundary Hodge and double-form operators.\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\\Riem{\\mathrm{Riem}}\n",
    "\\def\\sff{\\mathrm{I\\!I}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0941696",
   "metadata": {},
   "outputs": [],
   "source": [
    "import ngsolve\n",
    "from ngsolve import (\n",
    "    BBND,\n",
    "    BND,\n",
    "    VOL,\n",
    "    CF,\n",
    "    Cross,\n",
    "    Id,\n",
    "    InnerProduct,\n",
    "    Integrate,\n",
    "    Mesh,\n",
    "    Norm,\n",
    "    Normalize,\n",
    "    OuterProduct,\n",
    "    TaskManager,\n",
    "    specialcf,\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)\n",
    "\n",
    "\n",
    "def l2_error(a, b, vb=BND, bonus_intorder=4):\n",
    "    if vb == BND:\n",
    "        dX = ngsolve.dx(element_boundary=True, bonus_intorder=bonus_intorder)\n",
    "    elif vb == BBND:\n",
    "        dX = ngsolve.dx(element_vb=BBND, bonus_intorder=bonus_intorder)\n",
    "    else:\n",
    "        dX = ngsolve.dx(bonus_intorder=bonus_intorder)\n",
    "    return sqrt(Integrate(InnerProduct(a - b, a - b) * dX, mesh))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e519d4d",
   "metadata": {},
   "source": [
    "## Induced metric and projections\n",
    "\n",
    "Let $n$ be the $g$-unit normal and $n_E$ the Euclidean unit normal. The $g$-orthogonal and Euclidean tangential projectors are\n",
    "\n",
    "$$\n",
    "P_g=I-n\\otimes n^\\flat,\\qquad P_E=I-n_E\\otimes n_E.\n",
    "$$\n",
    "\n",
    "NGSDiffGeo embeds the induced metric and inverse metric into the ambient coordinates as\n",
    "\n",
    "$$\n",
    "g_F=g-n^\\flat\\otimes n^\\flat,\\qquad g_F^{-1}=g^{-1}-n\\otimes n.\n",
    "$$\n",
    "\n",
    "The corresponding properties are `G_F` and `G_F_inv`. On a codimension-two edge $E$, `G_E` and `G_E_inv` denote the metric and inverse metric induced along the $g$-unit edge tangent. Here the subscript `E` means *edge*, not Euclidean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84129cbb",
   "metadata": {},
   "outputs": [],
   "source": [
    "n = mf.normal\n",
    "n_E = specialcf.normal(3)\n",
    "P_g = Id(3) - OuterProduct(n, mf.Lower(n))\n",
    "P_E = Id(3) - OuterProduct(n_E, n_E)\n",
    "\n",
    "G_F = mf.G - OuterProduct(mf.Lower(n), mf.Lower(n))\n",
    "G_F_inv = mf.G_inv - OuterProduct(n, n)\n",
    "t = mf.tangent\n",
    "G_E = OuterProduct(mf.Lower(t), mf.Lower(t))\n",
    "G_E_inv = OuterProduct(t, t)\n",
    "\n",
    "assert l2_error(P_g, P_E * P_g) < TOL\n",
    "assert l2_error(P_E, P_g * P_E) < TOL\n",
    "assert l2_error(G_F, mf.G_F) < TOL\n",
    "assert l2_error(G_F_inv, mf.G_F_inv) < TOL\n",
    "assert l2_error(G_E, mf.G_E, vb=BBND) < TOL\n",
    "assert l2_error(G_E_inv, mf.G_E_inv, vb=BBND) < TOL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6106b79",
   "metadata": {},
   "source": [
    "## Codimension-two edge geometry\n",
    "\n",
    "On `BBND`, the two incident facets are indexed by $i=0,1$. The fields `edge_normal(i)` and `edge_conormal(i)` are their $g$-unit normals and their $g$-unit conormals within the facets. `edge_normals` and `edge_conormals` return both pairs. `ProjectTensor(..., \"E\")` and `ProjectDoubleForm(..., left=\"E\", right=\"E\")` project onto the edge tangent. A slot can instead be contracted with a chosen conormal using mode `\"m\"`.\n",
    "\n",
    "The metric line density on the edges is `VolumeForm(BBND)`. `AngleDefect` is the element-local Euclidean reference angle minus the angle measured with $g$. Both quantities are evaluated with `dx(element_vb=BBND)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f27ad36",
   "metadata": {},
   "outputs": [],
   "source": [
    "edge_normals = mf.edge_normals\n",
    "edge_conormals = mf.edge_conormals\n",
    "assert len(edge_normals) == len(edge_conormals) == 2\n",
    "\n",
    "for edge_normal, edge_conormal in zip(edge_normals, edge_conormals):\n",
    "    assert l2_error(mf.InnerProduct(t, edge_normal), CF(0), vb=BBND) < TOL\n",
    "    assert l2_error(mf.InnerProduct(t, edge_conormal), CF(0), vb=BBND) < TOL\n",
    "    assert l2_error(mf.InnerProduct(edge_normal, edge_conormal), CF(0), vb=BBND) < TOL\n",
    "    assert l2_error(mf.InnerProduct(edge_normal, edge_normal), CF(1), vb=BBND) < TOL\n",
    "    assert l2_error(mf.InnerProduct(edge_conormal, edge_conormal), CF(1), vb=BBND) < TOL\n",
    "\n",
    "v = dg.VectorField(CF((1 + x, y - 2, z + 3)))\n",
    "v_E = mf.ProjectTensor(v, \"E\")\n",
    "for edge_normal, edge_conormal in zip(edge_normals, edge_conormals):\n",
    "    assert l2_error(mf.InnerProduct(v_E, edge_normal), CF(0), vb=BBND) < TOL\n",
    "    assert l2_error(mf.InnerProduct(v_E, edge_conormal), CF(0), vb=BBND) < TOL\n",
    "\n",
    "alpha_E = mf.ProjectTensor(dg.OneForm(CF((1 + x, y, z**2))), \"E\")\n",
    "beta_E = mf.ProjectTensor(dg.OneForm(CF((z, 1 + y, x * y))), \"E\")\n",
    "phi_ambient = dg.DoubleForm(\n",
    "    dg.Einsum(\n",
    "        \"i,j->ij\", dg.OneForm(CF((1 + x, y, z**2))), dg.OneForm(CF((z, 1 + y, x * y)))\n",
    "    ),\n",
    "    p=1,\n",
    "    q=1,\n",
    "    dim=3,\n",
    ")\n",
    "phi_E = mf.ProjectDoubleForm(phi_ambient, left=\"E\", right=\"E\")\n",
    "phi_E_expected = dg.DoubleForm(dg.Einsum(\"i,j->ij\", alpha_E, beta_E), p=1, q=1, dim=3)\n",
    "assert l2_error(phi_E, phi_E_expected, vb=BBND) < TOL\n",
    "\n",
    "phi_m = mf.ProjectDoubleForm(\n",
    "    phi_ambient, left=\"m\", conormal=edge_conormals[0], project_remaining=False\n",
    ")\n",
    "assert (\n",
    "    l2_error(\n",
    "        phi_m, mf.ContractSlot(phi_ambient, edge_conormals[0], slot=\"left\"), vb=BBND\n",
    "    )\n",
    "    < TOL\n",
    ")\n",
    "\n",
    "dE = ngsolve.dx(element_vb=BBND, bonus_intorder=4)\n",
    "metric_edge_measure = Integrate(mf.VolumeForm(BBND) * dE, mesh)\n",
    "angle_defect_norm = sqrt(Integrate(mf.AngleDefect**2 * mf.VolumeForm(BBND) * dE, mesh))\n",
    "print(f\"metric edge measure: {metric_edge_measure:.6f}\")\n",
    "print(f\"weighted angle-defect norm: {angle_defect_norm:.6f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0be52bc",
   "metadata": {},
   "source": [
    "## Tangential connection and second fundamental form\n",
    "\n",
    "The boundary covariant derivative is selected with `vb=BND`. For a tangential vector field $X$, the Gauss formula is\n",
    "\n",
    "$$\n",
    "\\nabla_YX=\\nabla^F_YX+\\sff(Y,X)n.\n",
    "$$\n",
    "\n",
    "With NGSDiffGeo's convention, the second fundamental form is\n",
    "\n",
    "$$\n",
    "\\sff=-P_g^T(\\nabla n^\\flat)P_g.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffc973c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "sff_from_normal = -mf.ProjectTensor(mf.CovDerivative(mf.Lower(n)), \"F\")\n",
    "assert l2_error(sff_from_normal, mf.SFF) < TOL\n",
    "assert l2_error(mf.Trace(mf.SFF, vb=BND), mf.MeanCurvature) < TOL\n",
    "\n",
    "E1 = dg.VectorField(\n",
    "    Normalize(\n",
    "        ngsolve.IfPos(\n",
    "            Norm(n_E * CF((0, 0, 1))) - 0.9,\n",
    "            Cross(n_E, CF((1, 0, 0))),\n",
    "            Cross(n_E, CF((0, 0, 1))),\n",
    "        )\n",
    "    )\n",
    ")\n",
    "E2 = dg.VectorField(Cross(n_E, E1))\n",
    "X_ambient = (x - 0.1 * y**2 + 0.3 * x * z) * E1 + (y + 0.2 * z**2) * E2\n",
    "Y_ambient = (y - 0.2 * z**2) * E1 + (z + 0.1 * x**2) * E2\n",
    "X = mf.ProjectTensor(X_ambient, \"F\")\n",
    "Y = mf.ProjectTensor(Y_ambient, \"F\")\n",
    "\n",
    "ambient = mf.Contraction(mf.CovDerivative(X), Y)\n",
    "tangential_and_normal = (\n",
    "    mf.Contraction(mf.CovDerivative(X, vb=BND), Y) + mf.SFF[Y, X] * n\n",
    ")\n",
    "assert l2_error(ambient, tangential_and_normal) < TOL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99095d80",
   "metadata": {},
   "source": [
    "## Gauss and Codazzi equations\n",
    "\n",
    "The intrinsic curvature of $F$ is related to the ambient curvature by the Gauss equation\n",
    "\n",
    "$$\n",
    "\\Riem_F(X,Y,Z,W)=\\Riem(X,Y,Z,W)-\\tfrac12(\\sff\\owedge\\sff)(X,Y,Z,W).\n",
    "$$\n",
    "\n",
    "The mixed ambient curvature is described by the Codazzi equation. With the index and sign conventions used by NGSDiffGeo,\n",
    "\n",
    "$$\n",
    "\\Riem(X,Y,Z,n)=(\\nabla_X\\sff)(Y,Z)-(\\nabla_Y\\sff)(X,Z).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e67522c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "Z = (z + 0.3 * x * y - 0.1 * x**3) * E1 + (x - 0.1 * y**2) * E2\n",
    "\n",
    "curvature_term = mf.Riemann[X, Y, Z, n]\n",
    "cov_sff = mf.CovDerivative(mf.SFF)\n",
    "codazzi_term = cov_sff[X, Y, Z] - cov_sff[Y, X, Z]\n",
    "\n",
    "assert l2_error(curvature_term, codazzi_term) < TOL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c925a483",
   "metadata": {},
   "source": [
    "## Boundary Hodge star\n",
    "\n",
    "Passing `vb=BND` to `star` uses the induced metric on the two-dimensional boundary. For a tangential $k$-form $\\alpha$,\n",
    "\n",
    "$$\n",
    "\\star_F\\star_F\\alpha=(-1)^{k(2-k)}\\alpha,\\qquad\n",
    "\\alpha\\wedge\\star_F\\beta=\\langle\\alpha,\\beta\\rangle_{g_F}\\,\\omega_F.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86dfbc2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "one = dg.ScalarField(CF(1), dim=3)\n",
    "alpha = mf.ProjectTensor(dg.OneForm(CF((1 + x, y, z**2))), \"F\")\n",
    "beta = mf.ProjectTensor(dg.OneForm(CF((z, 1 + y, x * y))), \"F\")\n",
    "omega_F = mf.star(one, vb=BND)\n",
    "\n",
    "star_star_alpha = mf.star(mf.star(alpha, vb=BND), vb=BND)\n",
    "wedge_identity = dg.Wedge(alpha, mf.star(beta, vb=BND))\n",
    "inner_product_identity = mf.InnerProduct(alpha, beta, vb=BND, forms=True) * omega_F\n",
    "\n",
    "assert l2_error(star_star_alpha, -alpha) < TOL\n",
    "assert l2_error(wedge_identity, inner_product_identity) < TOL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1870fd92",
   "metadata": {},
   "source": [
    "## Exterior covariant derivatives on the boundary\n",
    "\n",
    "The same `vb=BND` convention selects the induced connection for double forms. Transposition exchanges the slots exactly as in the volume, and the coderivative can be recovered from the boundary Hodge star."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a80ff0f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "phi = dg.DoubleForm(dg.Einsum(\"i,j->ij\", alpha, beta), p=1, q=1, dim=3)\n",
    "phi = mf.ProjectDoubleForm(phi, left=\"F\", right=\"F\")\n",
    "\n",
    "d_right = mf.d_cov(phi, slot=\"right\", vb=BND)\n",
    "d_right_from_transpose = mf.d_cov(phi.trans, slot=\"left\", vb=BND).trans\n",
    "delta_right = mf.delta_cov(phi, slot=\"right\", vb=BND)\n",
    "delta_right_from_transpose = mf.delta_cov(phi.trans, slot=\"left\", vb=BND).trans\n",
    "\n",
    "assert l2_error(d_right, d_right_from_transpose) < TOL\n",
    "assert l2_error(delta_right, delta_right_from_transpose) < TOL\n",
    "\n",
    "surface_dim = 2\n",
    "sign = (-1) ** (surface_dim * phi.degree_left + surface_dim + 1)\n",
    "delta_from_star = sign * mf.star(\n",
    "    mf.d_cov(mf.star(phi, slot=\"left\", vb=BND), slot=\"left\", vb=BND),\n",
    "    slot=\"left\",\n",
    "    vb=BND,\n",
    ")\n",
    "assert l2_error(mf.delta_cov(phi, slot=\"left\", vb=BND), delta_from_star) < TOL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f71c8243",
   "metadata": {},
   "source": [
    "Finally, tangential projection does not commute with the ambient exterior covariant derivative. The correction is governed by the second fundamental form. For example,\n",
    "\n",
    "$$\n",
    "P^{tt}(d^\\nabla\\varphi)=d^{\\nabla^F}(P^{tt}\\varphi)-\\sff\\owedge P^{tn}\\varphi,\n",
    "$$\n",
    "\n",
    "with an analogous identity in the right slot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df4cb216",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha2 = dg.Wedge(\n",
    "    dg.OneForm(CF((x, 1 + y, z))),\n",
    "    dg.OneForm(CF((1 - z, x, 1 + y))),\n",
    ")\n",
    "phi12 = dg.DoubleForm(dg.Einsum(\"i,jk->ijk\", alpha, alpha2), p=1, q=2, dim=3)\n",
    "phi21 = dg.DoubleForm(dg.Einsum(\"ij,k->ijk\", alpha2, beta), p=2, q=1, dim=3)\n",
    "\n",
    "left = mf.ProjectDoubleForm(mf.d_cov(phi12, slot=\"left\"), left=\"F\", right=\"F\")\n",
    "right = mf.d_cov(\n",
    "    mf.ProjectDoubleForm(phi12, left=\"F\", right=\"F\"), slot=\"left\", vb=BND\n",
    ") - dg.Wedge(mf.SFF, mf.ProjectDoubleForm(phi12, left=\"F\", right=\"n\"))\n",
    "assert l2_error(left, right) < TOL\n",
    "\n",
    "left = mf.ProjectDoubleForm(mf.d_cov(phi21, slot=\"right\"), left=\"F\", right=\"F\")\n",
    "right = mf.d_cov(\n",
    "    mf.ProjectDoubleForm(phi21, left=\"F\", right=\"F\"), slot=\"right\", vb=BND\n",
    ") - dg.Wedge(mf.SFF, mf.ProjectDoubleForm(phi21, left=\"n\", right=\"F\"))\n",
    "assert l2_error(left, right) < TOL"
   ]
  }
 ],
 "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
}