From e4f3062a049c0436b18705d29ff3eb1ddfc6c0ad Mon Sep 17 00:00:00 2001 From: Luis Wirth Date: Thu, 23 Jul 2026 11:40:09 +0200 Subject: [PATCH] simplicial: the Levi-Civita connection of a Regge manifold Parallel transport across a facet is the unique isometry of the two cells' frames that restricts to the transition differential on the shared facet and carries the direction out of the source to the direction into the target: the unfolding of the pair into one flat frame. It exists because both cells read the facet's metric off the same edge lengths, so the connection is a derived quantity of the Regge primitive and needs no embedding. Piecewise flatness leaves no interior degrees of freedom and no contractible dual loop with holonomy, so the whole connection is one matrix per interior facet and all curvature concentrates on the codimension-2 hinges, where the holonomy rotates the normal plane by the deficit angle. This generalizes vertex_gaussian_curvature, whose 2D angle defect was the n = 2 case of deficit_angle; only the lumped-area density stays 2D, and the duplicate implementation is gone. The laws: metric compatibility, agreement of the shared facet metric, functoriality, trivial holonomy on a flat mesh, agreement of the holonomy rotation angle with the deficit angle, and Gauss-Bonnet. The sign of the unfolding is pinned against an embedding, where transport must be the ambient identity read in the two local frames -- a test, since the definition may not consult coordinates. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01DBMLrofgVMniugqfwti7p2 --- CLAUDE.md | 16 +- crates/simplicial/src/geometry.rs | 68 +- crates/simplicial/src/geometry/metric.rs | 2 + .../src/geometry/metric/connection.rs | 629 ++++++++++++++++++ crates/simplicial/src/topology/role.rs | 12 + 5 files changed, 689 insertions(+), 38 deletions(-) create mode 100644 crates/simplicial/src/geometry/metric/connection.rs diff --git a/CLAUDE.md b/CLAUDE.md index 5bbb8a5b..7f66913f 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -65,7 +65,7 @@ depending on nothing but `nalgebra-sparse`, joining the ladder only where | `gramian` | inner-product / metric structure | `Gramian` (non-degenerate symmetric, any signature), `Metric` (the pseudo-Riemannian metric tensor, any signature; Riemannian is $q = 0$), `CausalType` | | `coorder` | typed affine coordinates | `Coords` (coordinates tagged by their space), `affine::AffineTransform` | | `exterior` | the exterior algebra $Lambda^k$ | `ExteriorElement`, `Variance` (`Covariant`/`Contravariant`), `exterior_power`, wedge, interior product, musicals, Hodge star, `pullback`/`pushforward` of a value along a linear map | -| `simplicial` | the simplicial manifold $M_h$ | `topology::` (`Complex`, `Skeleton`, `SimplexRef`, the `role::` witnesses `Cell`/`Facet`/..., boundary operators, `orientation::Orientation`, `ordering::CellOrdering`, `refine::Subdivision`), `atlas::` (`Chart`, `MeshPoint`, `Transition`, `Bary`/`Local`, `SimplexQuadRule`), `geometry::` (`MeshLengthsSq` the intrinsic Regge primitive the engine consumes, `MeshCoords` and `CellGramians` the sources that convert into it) and `linalg::` (the dense/sparse nalgebra aliases and `CooMatrixExt` block-matrix builder every crate above it reuses) | +| `simplicial` | the simplicial manifold $M_h$ | `topology::` (`Complex`, `Skeleton`, `SimplexRef`, the `role::` witnesses `Cell`/`Facet`/..., boundary operators, `orientation::Orientation`, `ordering::CellOrdering`, `refine::Subdivision`), `atlas::` (`Chart`, `MeshPoint`, `Transition`, `Bary`/`Local`, `SimplexQuadRule`), `geometry::` (`MeshLengthsSq` the intrinsic Regge primitive the engine consumes, `MeshCoords` and `CellGramians` the sources that convert into it, `connection::Transport` the Levi-Civita connection and the hinge deficit angle) and `linalg::` (the dense/sparse nalgebra aliases and `CooMatrixExt` block-matrix builder every crate above it reuses) | | `glatt` | the continuum manifold $M$ | `Parametrization` (forward map $phi$, derived nearest-point chart, `sphere`/`ball`/`torus`/`graph`), `field::CoordField` (analytic data *on* $M$: `DiffFormClosure`, ...) | | `derham` | discrete differential forms | `Cochain`, `section::Section` (sections over the simplicial manifold) with the `Pullback` bridge (`pullback_on`/`pullback_through`) and `Sampler`, `interpolate::` (`WhitneyForm`, `WhitneyInterpolant`), `project::derham_map` | | `iterative` | matrix-free iterative solving | one object, an approximate inverse, reused as solver, preconditioner or smoother: stationary iteration, `Jacobi`, preconditioned `CG`, `MINRES` (symmetric indefinite), block-diagonal preconditioner; backend is `nalgebra-sparse` alone, no faer | @@ -191,6 +191,20 @@ and passes tests. pinning a metric accessor to the `Cell` witness would conflate *has a metric* (all simplices) with *carries a frame* (cells only). + This is also what makes the manifold have a **connection** at all. Two cells + sharing a facet read that facet's metric off the *same* edge lengths, so their + restrictions agree, and there is a unique isometry of their frames fixing the + facet and putting the two on opposite sides of it — the unfolding, + `MeshLengthsSq::transport`. A bag of unrelated `CellGramians` would admit no + such gluing, so the Levi-Civita connection is a *derived* quantity of the + Regge primitive and never an input. Being piecewise flat, the connection has + no degrees of freedom inside a cell and none across a contractible dual loop: + all of it lives on the facets, and all of its curvature on the codimension-2 + hinges, where the holonomy is a rotation by the deficit angle + (`Ridge::fan`, `deficit_angle`, `holonomy`). Curvature is a 2-form, hence + measured against area, which is *why* the hinges are codimension 2 — vertices + in 2D, edges in 3D, triangles in 4D, one mechanism and no special case. + A **point of the simplicial manifold** is therefore `MeshPoint` — a `Chart` plus barycentric coordinates — never a global coordinate, which on a Regge manifold does not exist. A **field** is a `Section`: a section of the diff --git a/crates/simplicial/src/geometry.rs b/crates/simplicial/src/geometry.rs index 76383d1a..b72a97d5 100644 --- a/crates/simplicial/src/geometry.rs +++ b/crates/simplicial/src/geometry.rs @@ -10,7 +10,11 @@ pub mod coord; pub mod metric; pub mod refine; -use crate::{Dim, atlas::refsimp_vol, topology::complex::Complex}; +use crate::{ + Dim, + atlas::refsimp_vol, + topology::{complex::Complex, role::roles}, +}; use gramian::Metric; @@ -27,25 +31,20 @@ pub fn cell_volume(metric: &Metric) -> f64 { } /// Discrete Gaussian curvature at every vertex of a 2-dimensional simplicial -/// manifold, by the angle defect: $K(v) = (2 pi - sum_(f ni v) theta_f (v)) \/ -/// A(v)$ at an interior vertex, or $(pi - sum_f theta_f (v)) \/ A(v)$ at a -/// boundary one -- the standard convention when a mesh has a rim, folding the -/// boundary's own geodesic curvature into $K$ rather than tracking it apart. -/// $A(v)$ is the barycentric lumped area $sum_(K ni v) "vol"(K) \/ 3$, the -/// standard mass-lumping convention. +/// manifold: the Regge deficit angle divided by the area it is spread over, +/// $K(v) = epsilon_v \/ A(v)$. /// -/// Intrinsic: reads the Regge edge lengths, not an embedding, since -/// [`Gramian::vertex_angle`](gramian::Gramian::vertex_angle) needs no -/// coordinates -- a Regge manifold given only as [`MeshLengthsSq`] has a -/// Gaussian curvature exactly as well as an embedded one, which is why the -/// primitive is what this consumes. This -/// Regge's curvature, concentrated at the codimension-2 hinges; in 2D the -/// hinges are vertices, which is the one case implemented here. Generalizing -/// to an $(n-2)$-dimensional hinge of an $n$-manifold needs dihedral angles -/// between codimension-1 facets, not corner angles between edges, and this -/// crate does not yet carry that computation -- fixed at 2D for the same -/// reason [`crate::mesher::sphere`] is: the concept itself, not a shortcut, is -/// what is 2-dimensional here. +/// Nothing here is 2-dimensional except the packaging. The deficit angle +/// ([`MeshLengthsSq::deficit_angle`]) is Regge curvature at a hinge in any +/// dimension, and a vertex is what a hinge is when $n = 2$; what this adds is +/// the density, $A(v) = sum_(K in v) vol(K) \/ 3$, the barycentric lumped area +/// under the standard mass-lumping convention. Gaussian curvature is a +/// *scalar field*, so it needs that division, and only in 2D is the deficit's +/// hinge a point for the density to sit at. +/// +/// Intrinsic: reads the Regge edge lengths, not an embedding -- a Regge +/// manifold given only as [`MeshLengthsSq`] has a Gaussian curvature exactly as +/// well as an embedded one, which is why the primitive is what this consumes. /// /// Exact, not an approximation of the smooth quantity: this is what /// Gauss-Bonnet defines discrete curvature to be, with @@ -58,30 +57,25 @@ pub fn vertex_gaussian_curvature(topology: &Complex, geometry: &MeshLengthsSq) - "Gaussian curvature is a 2D-surface quantity." ); let nvertices = topology.skeleton_raw(Dim::ZERO).len(); - let boundary: std::collections::HashSet = - topology.boundary_vertices().into_iter().collect(); - let mut angle_sum = vec![0.0; nvertices]; let mut areas = vec![0.0; nvertices]; for cell in topology.cells().handle_iter() { - let metric = geometry.cell_metric(cell); - let vol = cell_volume(&metric); - let verts = &cell.simplex().vertices; - for m in 0..3 { - let (a, b) = ((m + 1) % 3, (m + 2) % 3); - angle_sum[verts[m]] += metric.vector_gramian().vertex_angle(m, a, b); - areas[verts[m]] += vol / 3.0; + let vol = cell_volume(&geometry.cell_metric(cell)); + for &vertex in &cell.simplex().vertices { + areas[vertex] += vol / 3.0; } } - (0..nvertices) - .map(|v| { - let target = if boundary.contains(&v) { - std::f64::consts::PI - } else { - std::f64::consts::TAU - }; - (target - angle_sum[v]) / areas[v] + let hinges = topology + .role_skeleton::() + .expect("a 2-complex has ridges"); + hinges + .handle_iter() + .map(|hinge| { + let deficit = geometry + .deficit_angle(hinge) + .expect("a Riemannian surface has dihedral angles"); + deficit / areas[hinge.kidx()] }) .collect() } diff --git a/crates/simplicial/src/geometry/metric.rs b/crates/simplicial/src/geometry/metric.rs index c5377d58..15ad5ac9 100644 --- a/crates/simplicial/src/geometry/metric.rs +++ b/crates/simplicial/src/geometry/metric.rs @@ -1,7 +1,9 @@ +pub mod connection; pub mod geometry; pub mod mesh; pub mod simplex; +pub use connection::Transport; pub use geometry::CellGramians; pub type EdgeIdx = usize; diff --git a/crates/simplicial/src/geometry/metric/connection.rs b/crates/simplicial/src/geometry/metric/connection.rs new file mode 100644 index 00000000..15ce10ee --- /dev/null +++ b/crates/simplicial/src/geometry/metric/connection.rs @@ -0,0 +1,629 @@ +//! The Levi-Civita connection of a Regge manifold: parallel transport across +//! the facets, and the hinge holonomy that is its curvature. +//! +//! A connection is what identifies the tangent spaces at two points, which the +//! smooth structure alone does not: $T_p M$ and $T_q M$ are different vector +//! spaces with no canonical isomorphism between them. Given a metric there is +//! exactly one connection that is compatible with it (transport preserves the +//! inner product) and torsion-free -- the **Levi-Civita** connection -- and on +//! a Regge manifold it is not extra data at all, but a function of the squared +//! edge lengths. +//! +//! Piecewise flatness collapses the general machinery into finite linear +//! algebra, with no ODE and no Christoffel symbol anywhere: +//! +//! - **Inside a cell there is nothing to transport.** A cell is flat and +//! affine, so its chart's frame identifies all its tangent spaces at once, +//! exactly as in $RR^n$. The connection has no interior degrees of freedom. +//! - **All of the content sits on the facets.** Two cells sharing a facet +//! $sigma$ are two frames, and there is a unique isometry carrying one to the +//! other that fixes $sigma$ and puts the cells on opposite sides of it: the +//! *unfolding* of the pair into one flat frame. That is [`Transport`], and +//! the whole connection is one such map per interior facet. +//! - **Curvature lives on the hinges, and only there.** A loop in the dual +//! graph that avoids the codimension-2 skeleton bounds a region that unfolds +//! flat, so its transports telescope to the identity. What survives is the +//! holonomy around a *hinge* (a [`Ridge`]), an isometry fixing the hinge's +//! own tangent space and rotating the $2$-plane normal to it by the **deficit +//! angle** $epsilon_h = 2 pi - sum_(K supset h) theta_(K,h)$. +//! +//! Curvature is a $2$-form, so it is measured against area and never against a +//! length or a volume; a distribution supported on an $(n-2)$-simplex is +//! precisely a curvature $2$-form concentrated transversally to it. That is why +//! the hinges are the ridges, and it is why the deficit angle is the whole of +//! Regge curvature: vertices in 2D, edges in 3D, triangles in 4D, one mechanism. +//! +//! # Why the connection exists at all +//! +//! Two cells are not free to disagree arbitrarily on the facet they share: each +//! reads that facet's metric off the *same* squared edge lengths, so the two +//! restrictions agree ([`MeshLengthsSq::simplex_metric`] consults no containing +//! cell). Conformity of the geometry is what makes the gluing well posed, and a +//! bag of unrelated per-cell metrics would admit no connection at all. This is +//! the payoff of edge lengths being the primitive rather than +//! [`CellGramians`](super::CellGramians). +//! +//! # Frames, and what a transport matrix is +//! +//! Everything is expressed in each chart's own local cartesian frame -- the +//! basis $e_i = v_i - v_0$ of the cell's spanning vectors, in which +//! [`MeshLengthsSq::cell_metric`] is the Gramian. These frames are not +//! orthonormal, deliberately: orthonormalizing each cell would be an arbitrary +//! gauge choice per cell, whereas the local frame is canonical given the +//! chart. So a transport is not literally a matrix of $O(p,q)$ but an isometry +//! *between* two inner-product spaces, +//! +//! $ T^transpose g_(K') T = g_K, $ +//! +//! which is the same statement without a choice in it. A holonomy, whose source +//! and target are one chart, does land in the isometry group $O(g_K)$ of that +//! chart -- conjugate to $O(p,q)$, with the conjugation being exactly the gauge +//! that was not fixed. + +use super::mesh::MeshLengthsSq; +use crate::{ + Dim, + atlas::{Chart, ChartExt, ref_face_spanning_vectors, ref_vertices}, + topology::{ + handle::SimplexIdx, + role::{Cell, Ridge, roles}, + }, +}; + +use crate::linalg::{Matrix, Vector}; +use gramian::Metric; +use multiindex::Combination; + +/// Below this the $g$-norm of a would-be normal counts as null, and the +/// transverse direction fails to be a direction: the facet's induced metric is +/// degenerate and no isometry across it exists. +const NULL_EPS: f64 = 1e-12; + +/// Parallel transport between two charts: the linear isometry identifying the +/// frame of one cell with the frame of another. +/// +/// $ T_(K' arrow.l K): (RR^n, g_K) -> (RR^n, g_(K')), quad +/// T^transpose g_(K') T = g_K $ +/// +/// A vector of the manifold expressed in the source chart's frame, expressed +/// instead in the target's. Transports compose along a path +/// ([`then`](Self::then)) and invert along its reversal, so they are a functor +/// from the path groupoid of the dual graph into the isometries -- which is +/// the whole content of a connection. +#[derive(Debug, Clone)] +pub struct Transport { + source: SimplexIdx, + target: SimplexIdx, + matrix: Matrix, +} + +impl Transport { + /// The transport of a chart to itself: the identity, since a cell is flat + /// and its frame already identifies all of its tangent spaces. + pub fn identity(chart: Chart) -> Self { + let dim = chart.dim().index(); + Self { + source: chart.idx(), + target: chart.idx(), + matrix: Matrix::identity(dim, dim), + } + } + + pub fn source(&self) -> SimplexIdx { + self.source + } + pub fn target(&self) -> SimplexIdx { + self.target + } + pub fn dim(&self) -> Dim { + self.source.dim() + } + /// The matrix of the transport, in the local frames of the two charts. + pub fn matrix(&self) -> &Matrix { + &self.matrix + } + pub fn into_matrix(self) -> Matrix { + self.matrix + } + /// Whether source and target are the same chart. + pub fn is_identity(&self) -> bool { + self.source == self.target + } + + /// The same vector, in the target chart's frame. + pub fn apply(&self, vector: &Vector) -> Vector { + &self.matrix * vector + } + + /// The reverse transport, which is the inverse: transport along the reversed + /// path. An isometry is invertible, so this is total. + pub fn inverse(&self) -> Self { + Self { + source: self.target, + target: self.source, + matrix: self + .matrix + .clone() + .try_inverse() + .expect("a transport is an isometry, hence invertible"), + } + } + + /// Transport along the concatenation of two paths, $T_"next" compose T$. + /// + /// Panics if the paths do not meet: this one's target must be the next + /// one's source. + pub fn then(&self, next: &Self) -> Self { + assert_eq!( + self.target, next.source, + "Transports compose only along a connected path." + ); + Self { + source: self.source, + target: next.target, + matrix: &next.matrix * &self.matrix, + } + } +} + +/// The Levi-Civita connection, read off the Regge primitive. +impl MeshLengthsSq { + /// Parallel transport from one chart into an adjacent one, across the facet + /// they share. + /// + /// The unique linear isometry $(RR^n, g_K) -> (RR^n, g_(K'))$ that restricts + /// on the shared facet to the change of frame [`Transition::differential`] + /// and carries the direction *out of* the source to the direction *into* the + /// target -- that is, the unfolding of the two cells into one flat frame, + /// rather than the folding of one onto the other. The tangential condition + /// fixes it on an $(n-1)$-dimensional subspace, isometry fixes the remaining + /// direction up to sign, and the side condition picks the sign. + /// + /// Identity when source and target are the same chart. `None` when the two + /// are not adjacent, and `None` on the two degeneracies that make the + /// isometry not exist: a facet whose induced metric is degenerate (its normal + /// direction is null, so there is nothing to normalize), and two cells whose + /// metrics disagree in signature across the facet. + /// + /// [`Transition::differential`]: crate::atlas::Transition::differential + pub fn transport(&self, source: Chart, target: Chart) -> Option { + if source == target { + return Some(Transport::identity(source)); + } + // A 0-manifold has no facets, hence no adjacency; the total accessor is + // what says so, rather than a test on the dimension. + source.complex().role_skeleton::()?; + + let dim = source.dim(); + let shared = source.facets().find(|facet| { + let (a, b) = facet.adjacent_cells(); + a == target || b == Some(target) + })?; + + let source_positions = shared.simplex().relative_to(source.simplex()); + let target_positions = shared.simplex().relative_to(target.simplex()); + + // The tangent space of the facet, in each frame. The two are related by the + // transition differential, which *is* the change of frame there -- the + // metric-free half of the transport, already carried by the atlas. + let source_tangents = ref_face_spanning_vectors(dim, &source_positions); + let target_tangents = source.transition_to(target).differential() * &source_tangents; + + // Leaving the source through the facet is entering the target through it. + let source_metric = self.cell_metric(source); + let target_metric = self.cell_metric(target); + let leaving = outward_normal(&source_metric, dim, &source_positions)?; + let entering = -outward_normal(&target_metric, dim, &target_positions)?; + + // An isometry cannot exist if the transverse direction is spacelike on one + // side and timelike on the other. + if source_metric.norm_sq(&leaving).signum() != target_metric.norm_sq(&entering).signum() { + return None; + } + + let from = append_column(&source_tangents, &leaving); + let to = append_column(&target_tangents, &entering); + Some(Transport { + source: source.idx(), + target: target.idx(), + matrix: to * from.try_inverse()?, + }) + } + + /// Parallel transport along a path of charts, each adjacent to the next. + /// + /// The ordered product of the facet transports. `None` on an empty path, or + /// wherever a consecutive pair is not adjacent. + pub fn transport_along(&self, path: &[Chart]) -> Option { + let mut transport = Transport::identity(*path.first()?); + for step in path.windows(2) { + transport = transport.then(&self.transport(step[0], step[1])?); + } + Some(transport) + } + + /// The holonomy around a hinge: transport once around the fan of cells + /// meeting the ridge, back into the chart it started in. + /// + /// This is the curvature of the Regge manifold, in its integral form. It + /// fixes the hinge's own tangent space pointwise and acts on the $2$-plane + /// normal to it as a rotation by the deficit angle -- so it is the identity + /// exactly where the manifold is flat, and away from the hinges every loop + /// contracts and there is nothing else it could be. + /// + /// `None` on a boundary hinge, whose fan is open rather than closed and + /// around which there is no loop to transport along. + pub fn holonomy(&self, hinge: Ridge) -> Option { + if hinge.is_boundary() { + return None; + } + let fan = hinge.fan(); + let closed: Vec = fan.iter().copied().chain(std::iter::once(fan[0])).collect(); + self.transport_along(&closed) + } + + /// The rotation angle of the hinge holonomy, in $[0, pi]$. + /// + /// Read off the trace: the holonomy fixes the $(n-2)$-dimensional tangent + /// space of the hinge and rotates its normal plane, so + /// $tr H = (n - 2) + 2 cos epsilon$. This is the magnitude of the deficit + /// angle and not its sign, the sign being an orientation of the normal plane + /// that the trace cannot see; [`Self::deficit_angle`] carries the signed + /// value. + /// + /// Meaningful on a Riemannian geometry, where the normal plane is definite + /// and the holonomy is a genuine rotation. On an indefinite signature the + /// holonomy may be a boost, whose invariant is a rapidity rather than an + /// angle, and the returned value is then not it. + pub fn holonomy_angle(&self, hinge: Ridge) -> Option { + let holonomy = self.holonomy(hinge)?; + let dim = holonomy.dim().index() as f64; + let cos = 0.5 * (holonomy.matrix().trace() - (dim - 2.0)); + Some(cos.clamp(-1.0, 1.0).acos()) + } + + /// The interior dihedral angle a cell subtends at a hinge: the angle between + /// the two facets of the cell that contain the ridge. + /// + /// $ cos theta_(K,h) = - angle.l hat(n)_1, hat(n)_2 angle.r_g $ + /// + /// with $hat(n)_i$ the outward unit normals of those two facets, both read in + /// the cell's own frame from its Regge metric. In 2D, where a hinge is a + /// vertex, this is the corner angle of the triangle there. + /// + /// `None` when either facet's normal is null. Like + /// [`Self::holonomy_angle`], an angle is the Riemannian reading; on an + /// indefinite signature the quantity is a Lorentzian dihedral angle and the + /// arccosine is not it. + pub fn dihedral_angle(&self, cell: Cell, hinge: Ridge) -> Option { + let dim = cell.dim(); + let metric = self.cell_metric(cell); + let mut normals = cell + .facets() + .filter(|facet| hinge.simplex().is_subsimplex_of(facet.simplex())) + .map(|facet| { + let positions = facet.simplex().relative_to(cell.simplex()); + outward_normal(&metric, dim, &positions) + }); + let first = normals.next().expect("a cell has two facets at a hinge")?; + let second = normals.next().expect("a cell has two facets at a hinge")?; + Some((-metric.inner(&first, &second)).clamp(-1.0, 1.0).acos()) + } + + /// The deficit angle at a hinge: Regge curvature, as a scalar. + /// + /// $ epsilon_h = 2 pi - sum_(K supset h) theta_(K,h) $ + /// + /// The shortfall by which the cells around the hinge fail to close up flat -- + /// zero exactly when they do. On a boundary hinge the fan is open and the + /// closing target is $pi$ rather than $2 pi$, which folds the boundary's own + /// extrinsic bending into the same scalar rather than tracking it apart. + /// + /// Unlike [`Self::holonomy`] this needs no ordering of the fan: a sum is + /// commutative where a product of transports is not. It is signed, where + /// the holonomy's trace is not, and the two agree in magnitude. + pub fn deficit_angle(&self, hinge: Ridge) -> Option { + let target = if hinge.is_boundary() { + std::f64::consts::PI + } else { + std::f64::consts::TAU + }; + let sum: f64 = hinge + .get() + .cells() + .map(|cell| self.dihedral_angle(cell, hinge)) + .sum::>()?; + Some(target - sum) + } +} + +/// The $g$-unit normal of a facet within a cell, signed to point out of it. +/// +/// The normal direction is the $g$-orthogonal complement of the facet's tangent +/// space, which is one-dimensional and non-degenerate whenever the facet's +/// induced metric is; `None` is the null case, where there is no unit vector to +/// be had. Outwardness is the sign of the pairing with the vector to the +/// opposite vertex, which is non-zero precisely because that vector is +/// transverse. +fn outward_normal(metric: &Metric, dim: Dim, positions: &Combination) -> Option { + let tangents = ref_face_spanning_vectors(dim, positions); + let normal = unit_normal(metric, &tangents)?; + + let opposite = (0..=dim.index()) + .find(|&p| !positions.contains(p)) + .expect("a facet omits exactly one vertex of its cell"); + let vertices = ref_vertices(dim); + let inward = vertices.column(opposite) - vertices.column(positions.index_at(0)); + + Some(if metric.inner(&normal, &inward) > 0.0 { + -normal + } else { + normal + }) +} + +/// A vector spanning the $g$-orthogonal complement of the columns of +/// `tangents`, normalized to $g$-norm $plus.minus 1$. +/// +/// The complement is the kernel of $B^transpose g$, obtained as the eigenvector +/// of least eigenvalue of the normal equations -- which stays an $n times n$ +/// eigenproblem in every dimension, the empty tangent space of a $1$-manifold +/// included. `None` when the direction is null and admits no normalization. +fn unit_normal(metric: &Metric, tangents: &Matrix) -> Option { + let conditions = tangents.transpose() * metric.vector_gramian().matrix(); + let normal_equations = conditions.transpose() * &conditions; + + let eigen = na::SymmetricEigen::new(normal_equations); + let least = (0..eigen.eigenvalues.len()) + .min_by(|&a, &b| { + eigen.eigenvalues[a] + .abs() + .partial_cmp(&eigen.eigenvalues[b].abs()) + .unwrap() + }) + .expect("a cell of a manifold has at least one direction"); + let normal = eigen.eigenvectors.column(least).into_owned(); + + let norm_sq = metric.norm_sq(&normal); + (norm_sq.abs() > NULL_EPS).then(|| normal / norm_sq.abs().sqrt()) +} + +/// A matrix with one further column appended on the right. +fn append_column(matrix: &Matrix, column: &Vector) -> Matrix { + let mut extended = matrix.clone().insert_column(matrix.ncols(), 0.0); + extended.set_column(matrix.ncols(), column); + extended +} + +#[cfg(test)] +mod test { + use super::*; + use crate::{ + Dim, + geometry::coord::simplex::SimplexRefExt, + mesher::cartesian::CartesianGrid, + topology::{complex::Complex, role::roles}, + }; + + use approx::assert_relative_eq; + + fn adjacent_pairs(complex: &Complex) -> Vec<(Chart<'_>, Chart<'_>)> { + complex + .facets() + .handle_iter() + .filter_map(|facet| { + let (a, b) = facet.adjacent_cells(); + b.map(|b| (a, b)) + }) + .collect() + } + + /// Metric compatibility: $T^transpose g_(K') T = g_K$. Transport is an + /// isometry between the two frames, which is what makes it Levi-Civita + /// rather than merely a change of basis, and is the condition that pins the + /// transverse direction up to sign. + #[test] + fn transport_is_an_isometry() { + for dim in (1..=3usize).map(Dim::from) { + let (complex, coords) = CartesianGrid::new_unit(dim, 2).triangulate(); + let lengths = coords.to_edge_lengths_sq(&complex); + + for (source, target) in adjacent_pairs(&complex) { + let transport = lengths.transport(source, target).unwrap(); + let pullback = lengths + .cell_metric(target) + .vector_gramian() + .pullback(transport.matrix()); + assert_relative_eq!( + pullback.matrix(), + lengths.cell_metric(source).vector_gramian().matrix(), + epsilon = 1e-9 + ); + } + } + } + + /// The two cells induce the *same* metric on the facet they share -- they + /// read it off the same edge lengths -- which is why an isometry across the + /// facet exists at all, and why the connection is a function of the Regge + /// primitive and of nothing else. + #[test] + fn shared_facet_metric_agrees() { + for dim in (1..=3usize).map(Dim::from) { + let (complex, coords) = CartesianGrid::new_unit(dim, 2).triangulate(); + let lengths = coords.to_edge_lengths_sq(&complex); + + for facet in complex.facets().handle_iter() { + let (a, b) = facet.adjacent_cells(); + let Some(b) = b else { continue }; + let facet_metric = |cell: Chart| { + let positions = facet.simplex().relative_to(cell.simplex()); + let tangents = ref_face_spanning_vectors(cell.dim(), &positions); + lengths + .cell_metric(cell) + .vector_gramian() + .pullback(&tangents) + }; + assert_relative_eq!( + facet_metric(a).matrix(), + facet_metric(b).matrix(), + epsilon = 1e-12 + ); + } + } + } + + /// Against an embedding: on a mesh realized in $RR^n$ the unfolding of two + /// adjacent cells is the identity of the ambient space, so the transport in + /// the local frames must be $A_(K')^(-1) A_K$. + /// + /// This is the theorem that pins the construction, sign and all -- a + /// reflection would satisfy the isometry law just as well and fail here -- + /// and it is a *test* rather than the definition precisely because the + /// definition may not consult an embedding. + #[test] + fn transport_unfolds_an_embedded_mesh() { + for dim in (1..=3usize).map(Dim::from) { + let (complex, coords) = CartesianGrid::new_unit(dim, 2).triangulate(); + let lengths = coords.to_edge_lengths_sq(&complex); + + for (source, target) in adjacent_pairs(&complex) { + let source_frame = source.coord_simplex(&coords).linear_transform(); + let target_frame = target.coord_simplex(&coords).linear_transform(); + let expected = target_frame.try_inverse().unwrap() * source_frame; + + assert_relative_eq!( + lengths.transport(source, target).unwrap().matrix(), + &expected, + epsilon = 1e-9 + ); + } + } + } + + /// Transport reverses along the reversed path, and a chart transports to + /// itself as the identity: the path groupoid, discretely. + #[test] + fn transport_is_functorial() { + for dim in (1..=3usize).map(Dim::from) { + let (complex, coords) = CartesianGrid::new_unit(dim, 2).triangulate(); + let lengths = coords.to_edge_lengths_sq(&complex); + + for (source, target) in adjacent_pairs(&complex) { + let there = lengths.transport(source, target).unwrap(); + let back = lengths.transport(target, source).unwrap(); + assert_relative_eq!(back.matrix(), there.inverse().matrix(), epsilon = 1e-9); + + let roundtrip = there.then(&back); + assert!(roundtrip.is_identity()); + assert_relative_eq!( + roundtrip.matrix(), + &Matrix::identity(dim.index(), dim.index()), + epsilon = 1e-9 + ); + } + } + } + + /// A flat manifold has trivial holonomy around every hinge: the cells of a + /// triangulated box unfold into one frame, so every loop of the dual graph + /// telescopes. The base case of curvature, and the one every mesh with a + /// coordinate realization must pass. + #[test] + fn flat_mesh_has_no_holonomy() { + for dim in (2..=3usize).map(Dim::from) { + let (complex, coords) = CartesianGrid::new_unit(dim, 2).triangulate(); + let lengths = coords.to_edge_lengths_sq(&complex); + + let ridges = complex.role_skeleton::().unwrap(); + for hinge in ridges.handle_iter() { + let Some(holonomy) = lengths.holonomy(hinge) else { + continue; + }; + assert_relative_eq!( + holonomy.matrix(), + &Matrix::identity(dim.index(), dim.index()), + epsilon = 1e-8 + ); + assert_relative_eq!(lengths.holonomy_angle(hinge).unwrap(), 0.0, epsilon = 1e-6); + assert_relative_eq!(lengths.deficit_angle(hinge).unwrap(), 0.0, epsilon = 1e-9); + } + } + } + + /// The two readings of Regge curvature agree: the rotation angle of the + /// holonomy is the magnitude of the deficit angle. One is an ordered product + /// of isometries around the fan, the other a commutative sum of dihedral + /// angles, and they are the same number. + #[test] + fn holonomy_angle_is_the_deficit_angle() { + let (complex, coords) = crate::mesher::sphere::mesh_sphere_surface(2); + let lengths = coords.to_edge_lengths_sq(&complex); + + let ridges = complex.role_skeleton::().unwrap(); + let mut curved = 0; + for hinge in ridges.handle_iter() { + let Some(angle) = lengths.holonomy_angle(hinge) else { + continue; + }; + let deficit = lengths.deficit_angle(hinge).unwrap(); + assert_relative_eq!(angle, deficit.abs(), epsilon = 1e-8); + if deficit.abs() > 1e-6 { + curved += 1; + } + } + assert!(curved > 0, "a sphere is not flat"); + } + + /// Gauss-Bonnet on a closed surface: $sum_h epsilon_h = 2 pi chi$, the + /// deficit angles summing to the Euler characteristic with no area, no + /// refinement limit and no tolerance around a smooth quantity. On a sphere + /// $chi = 2$. + /// + /// The theorem that says the deficit angle *is* curvature rather than + /// merely resembling it. + #[test] + fn deficit_angles_sum_to_the_euler_characteristic() { + for refinement in 1..=3 { + let (complex, coords) = crate::mesher::sphere::mesh_sphere_surface(refinement); + let lengths = coords.to_edge_lengths_sq(&complex); + + let ridges = complex.role_skeleton::().unwrap(); + let total: f64 = ridges + .handle_iter() + .map(|hinge| lengths.deficit_angle(hinge).unwrap()) + .sum(); + + let euler: i64 = (0..=complex.dim().index()) + .map(|k| { + let n = complex.skeleton(k).len() as i64; + if k % 2 == 0 { n } else { -n } + }) + .sum(); + assert_relative_eq!(total, std::f64::consts::TAU * euler as f64, epsilon = 1e-9); + } + } + + /// In 2D a hinge is a vertex and a dihedral angle is a corner angle: the + /// general construction reproduces the elementary one it generalizes. + #[test] + fn dihedral_angle_in_2d_is_the_corner_angle() { + let (complex, coords) = CartesianGrid::new_unit(Dim::new(2), 2).triangulate(); + let lengths = coords.to_edge_lengths_sq(&complex); + + for cell in complex.cells().handle_iter() { + let metric = lengths.cell_metric(cell); + for (position, &vertex) in cell.simplex().vertices.iter().enumerate() { + let hinge = crate::topology::handle::SimplexIdx::new(Dim::ZERO, vertex) + .handle(&complex) + .role::(); + let (a, b) = ((position + 1) % 3, (position + 2) % 3); + assert_relative_eq!( + lengths.dihedral_angle(cell, hinge).unwrap(), + metric.vector_gramian().vertex_angle(position, a, b), + epsilon = 1e-9 + ); + } + } + } +} diff --git a/crates/simplicial/src/topology/role.rs b/crates/simplicial/src/topology/role.rs index 4bb0aa85..f9ccafc3 100644 --- a/crates/simplicial/src/topology/role.rs +++ b/crates/simplicial/src/topology/role.rs @@ -254,6 +254,18 @@ impl<'m> Ridge<'m> { .filter(move |facet| self.simplex().is_subsimplex_of(facet.simplex())) } + /// Whether the hinge lies on the boundary: one of the facets containing it + /// bounds a single cell, so its [`fan`](Self::fan) is open rather than + /// closed. The codimension-2 reading of [`Facet::is_boundary`], and the + /// condition under which there is no loop around the hinge to transport + /// along. + pub fn is_boundary(self) -> bool { + self + .get() + .cells() + .any(|cell| self.hinge_facets(cell).any(Facet::is_boundary)) + } + /// The fan of the hinge: the incident cells in adjacency order around the /// ridge, consecutive ones sharing a facet that contains it. Closed (the /// last cell neighboring the first) iff the ridge is interior; an open fan