Differential operators and curvature on hypersurfaces

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.

\[\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}}} \def\Riem{\mathrm{Riem}} \def\sff{\mathrm{I\!I}}\]

[1]:

import ngsolve
from ngsolve import (
    BBND,
    BND,
    VOL,
    CF,
    Cross,
    Id,
    InnerProduct,
    Integrate,
    Mesh,
    Norm,
    Normalize,
    OuterProduct,
    TaskManager,
    specialcf,
    sqrt,
    x,
    y,
    z,
)
from netgen.occ import unit_cube
import ngsdiffgeo as dg

TOL = 1e-7
mesh = Mesh(unit_cube.GenerateMesh(maxh=2))
mf = dg.RiemannianManifold(dg.Heisenberg().metric)


def l2_error(a, b, vb=BND, bonus_intorder=4):
    if vb == BND:
        dX = ngsolve.dx(element_boundary=True, bonus_intorder=bonus_intorder)
    elif vb == BBND:
        dX = ngsolve.dx(element_vb=BBND, bonus_intorder=bonus_intorder)
    else:
        dX = ngsolve.dx(bonus_intorder=bonus_intorder)
    return sqrt(Integrate(InnerProduct(a - b, a - b) * dX, mesh))

Induced metric and projections

Let \(n\) be the \(g\)-unit normal and \(n_E\) the Euclidean unit normal. The \(g\)-orthogonal and Euclidean tangential projectors are

\[P_g=I-n\otimes n^\flat,\qquad P_E=I-n_E\otimes n_E.\]

NGSDiffGeo embeds the induced metric and inverse metric into the ambient coordinates as

\[g_F=g-n^\flat\otimes n^\flat,\qquad g_F^{-1}=g^{-1}-n\otimes 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.

[2]:

n = mf.normal
n_E = specialcf.normal(3)
P_g = Id(3) - OuterProduct(n, mf.Lower(n))
P_E = Id(3) - OuterProduct(n_E, n_E)

G_F = mf.G - OuterProduct(mf.Lower(n), mf.Lower(n))
G_F_inv = mf.G_inv - OuterProduct(n, n)
t = mf.tangent
G_E = OuterProduct(mf.Lower(t), mf.Lower(t))
G_E_inv = OuterProduct(t, t)

assert l2_error(P_g, P_E * P_g) < TOL
assert l2_error(P_E, P_g * P_E) < TOL
assert l2_error(G_F, mf.G_F) < TOL
assert l2_error(G_F_inv, mf.G_F_inv) < TOL
assert l2_error(G_E, mf.G_E, vb=BBND) < TOL
assert l2_error(G_E_inv, mf.G_E_inv, vb=BBND) < TOL

Codimension-two edge geometry

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".

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).

[3]:

edge_normals = mf.edge_normals
edge_conormals = mf.edge_conormals
assert len(edge_normals) == len(edge_conormals) == 2

for edge_normal, edge_conormal in zip(edge_normals, edge_conormals):
    assert l2_error(mf.InnerProduct(t, edge_normal), CF(0), vb=BBND) < TOL
    assert l2_error(mf.InnerProduct(t, edge_conormal), CF(0), vb=BBND) < TOL
    assert l2_error(mf.InnerProduct(edge_normal, edge_conormal), CF(0), vb=BBND) < TOL
    assert l2_error(mf.InnerProduct(edge_normal, edge_normal), CF(1), vb=BBND) < TOL
    assert l2_error(mf.InnerProduct(edge_conormal, edge_conormal), CF(1), vb=BBND) < TOL

v = dg.VectorField(CF((1 + x, y - 2, z + 3)))
v_E = mf.ProjectTensor(v, "E")
for edge_normal, edge_conormal in zip(edge_normals, edge_conormals):
    assert l2_error(mf.InnerProduct(v_E, edge_normal), CF(0), vb=BBND) < TOL
    assert l2_error(mf.InnerProduct(v_E, edge_conormal), CF(0), vb=BBND) < TOL

alpha_E = mf.ProjectTensor(dg.OneForm(CF((1 + x, y, z**2))), "E")
beta_E = mf.ProjectTensor(dg.OneForm(CF((z, 1 + y, x * y))), "E")
phi_ambient = dg.DoubleForm(
    dg.Einsum(
        "i,j->ij", dg.OneForm(CF((1 + x, y, z**2))), dg.OneForm(CF((z, 1 + y, x * y)))
    ),
    p=1,
    q=1,
    dim=3,
)
phi_E = mf.ProjectDoubleForm(phi_ambient, left="E", right="E")
phi_E_expected = dg.DoubleForm(dg.Einsum("i,j->ij", alpha_E, beta_E), p=1, q=1, dim=3)
assert l2_error(phi_E, phi_E_expected, vb=BBND) < TOL

phi_m = mf.ProjectDoubleForm(
    phi_ambient, left="m", conormal=edge_conormals[0], project_remaining=False
)
assert (
    l2_error(
        phi_m, mf.ContractSlot(phi_ambient, edge_conormals[0], slot="left"), vb=BBND
    )
    < TOL
)

dE = ngsolve.dx(element_vb=BBND, bonus_intorder=4)
metric_edge_measure = Integrate(mf.VolumeForm(BBND) * dE, mesh)
angle_defect_norm = sqrt(Integrate(mf.AngleDefect**2 * mf.VolumeForm(BBND) * dE, mesh))
print(f"metric edge measure: {metric_edge_measure:.6f}")
print(f"weighted angle-defect norm: {angle_defect_norm:.6f}")

metric edge measure: 77.335631
weighted angle-defect norm: 1.904668

Tangential connection and second fundamental form

The boundary covariant derivative is selected with vb=BND. For a tangential vector field \(X\), the Gauss formula is

\[\nabla_YX=\nabla^F_YX+\sff(Y,X)n.\]

With NGSDiffGeo’s convention, the second fundamental form is

\[\sff=-P_g^T(\nabla n^\flat)P_g.\]

[4]:

sff_from_normal = -mf.ProjectTensor(mf.CovDerivative(mf.Lower(n)), "F")
assert l2_error(sff_from_normal, mf.SFF) < TOL
assert l2_error(mf.Trace(mf.SFF, vb=BND), mf.MeanCurvature) < TOL

E1 = dg.VectorField(
    Normalize(
        ngsolve.IfPos(
            Norm(n_E * CF((0, 0, 1))) - 0.9,
            Cross(n_E, CF((1, 0, 0))),
            Cross(n_E, CF((0, 0, 1))),
        )
    )
)
E2 = dg.VectorField(Cross(n_E, E1))
X_ambient = (x - 0.1 * y**2 + 0.3 * x * z) * E1 + (y + 0.2 * z**2) * E2
Y_ambient = (y - 0.2 * z**2) * E1 + (z + 0.1 * x**2) * E2
X = mf.ProjectTensor(X_ambient, "F")
Y = mf.ProjectTensor(Y_ambient, "F")

ambient = mf.Contraction(mf.CovDerivative(X), Y)
tangential_and_normal = (
    mf.Contraction(mf.CovDerivative(X, vb=BND), Y) + mf.SFF[Y, X] * n
)
assert l2_error(ambient, tangential_and_normal) < TOL

Gauss and Codazzi equations

The intrinsic curvature of \(F\) is related to the ambient curvature by the Gauss equation

\[\Riem_F(X,Y,Z,W)=\Riem(X,Y,Z,W)-\tfrac12(\sff\owedge\sff)(X,Y,Z,W).\]

The mixed ambient curvature is described by the Codazzi equation. With the index and sign conventions used by NGSDiffGeo,

\[\Riem(X,Y,Z,n)=(\nabla_X\sff)(Y,Z)-(\nabla_Y\sff)(X,Z).\]

[5]:

Z = (z + 0.3 * x * y - 0.1 * x**3) * E1 + (x - 0.1 * y**2) * E2

curvature_term = mf.Riemann[X, Y, Z, n]
cov_sff = mf.CovDerivative(mf.SFF)
codazzi_term = cov_sff[X, Y, Z] - cov_sff[Y, X, Z]

assert l2_error(curvature_term, codazzi_term) < TOL

Boundary Hodge star

Passing vb=BND to star uses the induced metric on the two-dimensional boundary. For a tangential \(k\)-form \(\alpha\),

\[\star_F\star_F\alpha=(-1)^{k(2-k)}\alpha,\qquad \alpha\wedge\star_F\beta=\langle\alpha,\beta\rangle_{g_F}\,\omega_F.\]

[6]:

one = dg.ScalarField(CF(1), dim=3)
alpha = mf.ProjectTensor(dg.OneForm(CF((1 + x, y, z**2))), "F")
beta = mf.ProjectTensor(dg.OneForm(CF((z, 1 + y, x * y))), "F")
omega_F = mf.star(one, vb=BND)

star_star_alpha = mf.star(mf.star(alpha, vb=BND), vb=BND)
wedge_identity = dg.Wedge(alpha, mf.star(beta, vb=BND))
inner_product_identity = mf.InnerProduct(alpha, beta, vb=BND, forms=True) * omega_F

assert l2_error(star_star_alpha, -alpha) < TOL
assert l2_error(wedge_identity, inner_product_identity) < TOL

Exterior covariant derivatives on the boundary

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.

[7]:

phi = dg.DoubleForm(dg.Einsum("i,j->ij", alpha, beta), p=1, q=1, dim=3)
phi = mf.ProjectDoubleForm(phi, left="F", right="F")

d_right = mf.d_cov(phi, slot="right", vb=BND)
d_right_from_transpose = mf.d_cov(phi.trans, slot="left", vb=BND).trans
delta_right = mf.delta_cov(phi, slot="right", vb=BND)
delta_right_from_transpose = mf.delta_cov(phi.trans, slot="left", vb=BND).trans

assert l2_error(d_right, d_right_from_transpose) < TOL
assert l2_error(delta_right, delta_right_from_transpose) < TOL

surface_dim = 2
sign = (-1) ** (surface_dim * phi.degree_left + surface_dim + 1)
delta_from_star = sign * mf.star(
    mf.d_cov(mf.star(phi, slot="left", vb=BND), slot="left", vb=BND),
    slot="left",
    vb=BND,
)
assert l2_error(mf.delta_cov(phi, slot="left", vb=BND), delta_from_star) < TOL

Finally, tangential projection does not commute with the ambient exterior covariant derivative. The correction is governed by the second fundamental form. For example,

\[P^{tt}(d^\nabla\varphi)=d^{\nabla^F}(P^{tt}\varphi)-\sff\owedge P^{tn}\varphi,\]

with an analogous identity in the right slot.

[8]:

alpha2 = dg.Wedge(
    dg.OneForm(CF((x, 1 + y, z))),
    dg.OneForm(CF((1 - z, x, 1 + y))),
)
phi12 = dg.DoubleForm(dg.Einsum("i,jk->ijk", alpha, alpha2), p=1, q=2, dim=3)
phi21 = dg.DoubleForm(dg.Einsum("ij,k->ijk", alpha2, beta), p=2, q=1, dim=3)

left = mf.ProjectDoubleForm(mf.d_cov(phi12, slot="left"), left="F", right="F")
right = mf.d_cov(
    mf.ProjectDoubleForm(phi12, left="F", right="F"), slot="left", vb=BND
) - dg.Wedge(mf.SFF, mf.ProjectDoubleForm(phi12, left="F", right="n"))
assert l2_error(left, right) < TOL

left = mf.ProjectDoubleForm(mf.d_cov(phi21, slot="right"), left="F", right="F")
right = mf.d_cov(
    mf.ProjectDoubleForm(phi21, left="F", right="F"), slot="right", vb=BND
) - dg.Wedge(mf.SFF, mf.ProjectDoubleForm(phi21, left="n", right="F"))
assert l2_error(left, right) < TOL