∑MathAdvanced

Curvature

Key Points

•
Curvature measures how a geometric object bends, and it comes in several flavors: Gaussian, sectional, and Ricci curvature.
•
Gaussian curvature captures intrinsic bending of a surface at a point using the product of principal curvatures, independent of how the surface sits in space.
•
Sectional curvature measures the curvature of a two-dimensional plane slice through a higher-dimensional manifold.
•
Ricci curvature aggregates sectional curvatures around a direction to describe average volume distortion and geodesic focusing.
•
On a sphere of radius R, all sectional curvatures are 1/ $R^{2}$ , while on a hyperbolic plane they are -1/ $R^{2}$ .
•
Gaussian curvature can be computed from first and second fundamental forms or, for z = f(x, y), from derivatives of f using a closed-form formula.
•
In 2D, the Ricci scalar equals 2 times the Gaussian curvature, making computations notably simpler.
•
Discrete meshes approximate curvature using angle deficits at vertices and local area estimates, enabling practical computation in C++.

Prerequisites

→Multivariable calculus — Curvature formulas require partial derivatives, gradients, Hessians, and Taylor approximations.
→Linear algebra — Metrics, Christoffel symbols, and curvature tensors rely on matrix operations and inner products.
→Vector calculus and differential geometry basics — Understanding tangent spaces, normal vectors, and geodesics is essential for definitions of curvature.
→Numerical analysis — Stable finite differences and error control are needed for accurate curvature computation.
→C++ programming fundamentals — Implementations require arrays, functions, and careful floating-point handling.

Detailed Explanation

Tap terms for definitions

01Overview

Curvature generalizes the idea of how a curve bends to higher-dimensional shapes like surfaces and abstract spaces called manifolds. For surfaces in three-dimensional space, the Gaussian curvature at a point reflects how the surface bends in two principal (mutually orthogonal) directions; it is positive on dome-like regions (e.g., a sphere), negative on saddle-like regions, and zero on flat regions. In higher dimensions, sectional curvature measures curvature on a chosen two-dimensional plane within the tangent space, allowing us to probe the manifold’s bending in various directions. Ricci curvature condenses this directional information into an average over planes containing a given direction; it is central in geometry and physics because it influences how volumes change along geodesics and features in Einstein’s field equations. Curvature can be computed from local differential data: the metric tensor (which encodes lengths and angles), its derivatives (via Christoffel symbols), and combinations that form the Riemann curvature tensor. For surfaces given as graphs z = f(x, y), formulas using first and second derivatives of f give Gaussian curvature directly. On triangle meshes, discrete approximations use angle deficits and local area to estimate curvature numerically. This topic bridges pure mathematics, computer graphics, simulation, and general relativity, and can be explored effectively with C++ implementations using numerical differentiation and basic linear algebra.

02Intuition & Analogies

Imagine driving over different landscapes. On a perfectly flat parking lot, you never feel the car tilt in any direction—this is zero curvature. On a spherical hilltop (like the top of a dome), no matter which way you face, the ground curves downward; that omnidirectional bending is positive curvature. On a saddle-shaped surface, if you face one way the ground curves up, but turn 90 degrees and it curves down—that mixed behavior signals negative curvature. Gaussian curvature captures this by multiplying the maximum and minimum bendings (principal curvatures); same sign means dome-like or bowl-like, opposite signs mean saddle. For higher dimensions, think of slicing an object with a thin sheet of paper. Each slice gives a 2D cross-section. The curvature you see in that slice is the sectional curvature for the plane of the slice. By checking many slices, you understand the object’s overall bending behavior. Ricci curvature is like asking: if I release a cloud of tiny particles moving straight ahead (along geodesics), does the cloud tend to focus or spread? Positive Ricci curvature nudges geodesics together (focusing, like gravity near a massive body), while negative Ricci curvature lets them drift apart. In two dimensions, this average behavior is entirely determined by the Gaussian curvature, but in higher dimensions, Ricci summarizes many sectional curvatures at once. In computation, we approximate these ideas with derivatives and linear algebra. For a height field z = f(x, y), curvature comes from how quickly slopes change. For arbitrary coordinates, we use the metric tensor to measure lengths/angles and compute connection coefficients (Christoffel symbols), which describe how vectors change as we move. Combining their derivatives produces the Riemann tensor, from which sectional and Ricci curvatures are extracted.

03Formal Definition

Let (M, g) be a smooth n-dimensional Riemannian manifold with metric tensor g. The Levi–Civita connection ∇ is the unique torsion-free connection compatible with g. The Riemann curvature tensor is the (1,3)-tensor defined by R(X, Y)Z = ∇_X ∇_Y Z − ∇_Y ∇_X Z − ∇_[X,Y] Z. In coordinates, its components are

R_{ijk}^{ℓ}

\partial_{j}

Γ_{ik}^{ℓ}

−

\partial_{i}

Γ_{jk}^{ℓ}

Γ_{jm}^{ℓ}

Γ_{ik}^{m}

−

Γ_{im}^{ℓ}

Γ_{jk}^{m}

, where

Γ_{ij}^{k}

are Christoffel symbols and indices repeat under summation. For a 2-plane

σ

span

\{u, v\}

\subset

T_{p}

M, the sectional curvature is K(

σ

) =

\frac{⟨ R ( u , v ) v , u ⟩}{∥ u ∥ ^{2} ∥ v ∥ ^{2} - ⟨ u , v ⟩ ^{2}}

. The Ricci curvature is the trace of R on the first and third indices:

Ric

(X, Y) =

tr

\mapsto

R(Z, X)Y), or in coordinates

Ric_{ij}

R_{ikj}^{k}

. The scalar curvature is

R = g^{ij}

Ric_{ij}

. For an embedded regular surface S

\subset

R^{3}

with first fundamental form I and second fundamental form II, the Gaussian curvature is K =

\frac{d e t ( II )}{d e t ( I )}

k_{1}

k_{2}

, where

k_{1}

and

k_{2}

are principal curvatures. For a Monge patch

z = f (x

, y), K =

\frac{f _{xx} f _{yy} - f _{x y}^{2}}{( 1 + f _{x}^{2} + f _{y}^{2} ) ^{2}}

. These definitions are intrinsic (depend only on g) except for expressions using embeddings, which are equal to the intrinsic ones by theorems of Gauss.

04When to Use

Use Gaussian curvature when analyzing 2D surfaces, such as terrain models, CAD surfaces, and mesh processing. It distinguishes spherical-like bulges (K > 0), saddle regions (K < 0), and flat areas (K ≈ 0), guiding smoothing, feature detection, and remeshing. Use sectional curvature in higher-dimensional geometry when you need direction-dependent curvature on a manifold. It is essential in understanding geodesic behavior, comparison theorems, and the geometry of spaces with constant curvature (spherical, Euclidean, hyperbolic). In algorithms, it appears when you model data on manifolds (e.g., optimization on matrix manifolds) and want to reason about local geometry. Use Ricci curvature when you care about average focusing/defocusing along directions or volume growth—e.g., in general relativity (Einstein equations relate Ricci to matter), heat kernel estimates, and geometric analysis. In data science, discrete analogs (Ollivier/Ricci flow on graphs) help detect community structure and robustness. In 2D, Ricci flow uniformizes metrics; in graphics, discrete Ricci flow aids mesh parameterization. Computationally: choose closed-form Gaussian curvature from derivatives for height fields; choose discrete angle-deficit curvature for triangle meshes; choose Christoffel/Riemann-based computations for general coordinates or when you have an explicit metric (e.g., sphere in spherical coordinates).

⚠️Common Mistakes

Mixing intrinsic and extrinsic viewpoints: Gaussian curvature is intrinsic; it does not change if you bend a surface without stretching (e.g., rolling paper into a cylinder). Avoid formulas that depend on embedding when you only have a metric.
Numerical differentiation pitfalls: Finite difference step sizes that are too large (biased) or too small (catastrophic cancellation) can destroy accuracy in curvature formulas involving second derivatives. Use central differences and tune step sizes; consider automatic differentiation when possible.
Coordinate singularities: Spherical coordinates are singular at the poles; naive computations of Christoffel symbols and curvature can blow up numerically even though true curvature is finite. Work away from singularities or change charts.
Index mistakes in tensor formulas: Confusing positions of indices (upper/lower) or summation ranges leads to sign or factor errors. Write out small cases (e.g., 2D) and verify symmetries R_{ijkl} = −R_{jikl} = −R_{ijlk} = R_{klij}.
Using non-orthonormal bases in sectional curvature without the denominator: The formula requires normalization by |u|^2 |v|^2 − \langle u, v \rangle^2; forgetting this yields scale-dependent results.
Mesh curvature without proper area: Angle deficit divided by an arbitrary area gives meaningless magnitudes. Use a consistent area (e.g., 1/3 of incident triangle areas or mixed Voronoi area) and ensure triangles are well-shaped to reduce discretization error.

Key Formulas

Gaussian Curvature via Principal Curvatures

K = k_{1} k_{2}

Explanation: The Gaussian curvature equals the product of the maximum and minimum normal curvatures at a point. Positive means dome-like, negative means saddle-like, and zero means flat in at least one principal direction.

Gauss Formula (First/Second Fundamental Forms)

K = \frac{e g - f ^{2}}{EG - F ^{2}}

Explanation: For a parametric surface with first fundamental form entries E, F, G and second fundamental form entries e, f, g, the Gaussian curvature is the ratio of determinants. It is intrinsic even though each term is defined extrinsically.

Monge Patch Gaussian Curvature

K = \frac{f _{xx} f _{yy} - f _{x y}^{2}}{( 1 + f _{x}^{2} + f _{y}^{2} ) ^{2}}

Explanation: For a surface z = f(x, y), Gaussian curvature depends on the Hessian of f and the gradient magnitude. It is convenient for height fields where derivatives can be computed numerically.

Christoffel Symbols

Γ_{ij}^{k} = \frac{1}{2} g^{k ℓ} (\partial_{i} g_{ℓ j} + \partial_{j} g_{ℓ i} - \partial_{ℓ} g_{ij})

Explanation: Connection coefficients computed from the metric tensor and its first derivatives. They describe how basis vectors change from point to point and enter geodesic and curvature computations.

Riemann Curvature Tensor (Components)

R_{ijk}^{ℓ} = \partial_{j} Γ_{ik}^{ℓ} - \partial_{i} Γ_{jk}^{ℓ} + Γ_{jm}^{ℓ} Γ_{ik}^{m} - Γ_{im}^{ℓ} Γ_{jk}^{m}

Explanation: This combines derivatives of Christoffel symbols and quadratic terms in them. It measures the intrinsic curvature independent of coordinates.

Sectional Curvature

K (u, v) = \frac{⟨ R ( u , v ) v , u ⟩}{∥ u ∥ ^{2} ∥ v ∥ ^{2} - ⟨ u , v ⟩ ^{2}}

Explanation: Curvature of the 2D plane spanned by u and v. The denominator normalizes for the area of the parallelogram formed by u and v.

Ricci Curvature (Contraction)

Ric_{ij} = R_{ikj}^{k}

Explanation: Ricci curvature is obtained by tracing the Riemann tensor on the first and third indices. It summarizes average sectional curvatures through each direction.

Scalar Curvature

R = g^{ij} Ric_{ij}

Explanation: The scalar curvature is the trace of the Ricci tensor with the inverse metric. In two dimensions, R = 2K.

2D Curvature from Riemann

K = \frac{R _{1212}}{det ( g )} (n = 2)

Explanation: In two dimensions, there is only one independent sectional curvature and it equals the Gaussian curvature. It can be computed from the Riemann tensor and the determinant of the metric.

Gauss–Bonnet Theorem

\int_{M} K d A = 2 π χ (M)

Explanation: The total Gaussian curvature of a compact 2D surface equals 2π times its Euler characteristic. This links geometry to topology.

Constant Curvature Models

K_{S^{n}} = \frac{1}{R ^{2}}, K_{H^{n}} = - \frac{1}{R ^{2}}

Explanation: Spheres have constant positive sectional curvature 1/ $R^{2}$ , while hyperbolic spaces have constant negative curvature −1/ $R^{2}$ . These are canonical examples for intuition and testing.

Complexity Analysis

Curvature computations vary in cost depending on representation and discretization. - Monge patch (

z = f (x

, y)): At a single point, evaluating Gaussian curvature using central finite differences requires O(1) calls to f to estimate first and second derivatives (typically 6–9 evaluations). The arithmetic is constant-time, so time is O(1) per point and O(NM) for an N×M grid. Space is O(1) per point; when storing a full curvature grid, space is O(NM). - Metric-based (2D) curvature via Christoffel/Riemann: Computing the metric, its inverse and determinant is O(1) in 2D. Estimating Christoffel symbols requires O(1) metric evaluations and arithmetic. To estimate Riemann components, we need derivatives of Christoffel symbols, which we approximate by recomputing Christoffels at nearby points; this still stays O(1) per evaluation. Numerical stability depends on the step size; too small increases rounding error, too large increases truncation error. Space is O(1). - Discrete mesh angle-deficit: For a mesh with V vertices, E edges, and F faces, accumulating angles per vertex visits each face once and does O(1) work per corner (three corners per triangle). Time is O(F) and often O(V) for well-shaped meshes (since F ≈ O(V)). Memory is O(V + F) to store the mesh and O(V) for curvature values. The algorithm’s accuracy improves with mesh refinement and well-shaped (non-sliver) triangles. In all cases, the dominant costs are evaluations of user-provided functions (e.g., f in Monge patch) and basic linear algebra (2×2 matrix inverses, dot/cross products). For higher-dimensional curvature (

n > 2

), naive tensor computations scale as O(

n^{4}

) to assemble all Riemann components, but many symmetries reduce actual cost; this answer focuses on 2D where costs are constant per point.

Code Examples

Gaussian curvature of a height field z = f(x, y) via central differences

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 // Example height function: a saddle z = (x^2 - y^2)/2
5 static double f(double x, double y) {
6     return 0.5 * (x*x - y*y);
7 }
8 
9 // Central finite differences for first and second derivatives
10 struct Derivs {
11     double fx, fy, fxx, fyy, fxy;
12 };
13 
14 Derivs derivatives(double x, double y, double h = 1e-4) {
15     Derivs d{};
16     // First derivatives
17     d.fx = (f(x + h, y) - f(x - h, y)) / (2*h);
18     d.fy = (f(x, y + h) - f(x, y - h)) / (2*h);
19     // Second derivatives
20     d.fxx = (f(x + h, y) - 2*f(x, y) + f(x - h, y)) / (h*h);
21     d.fyy = (f(x, y + h) - 2*f(x, y) + f(x, y - h)) / (h*h);
22     // Mixed derivative via central differences
23     d.fxy = (f(x + h, y + h) - f(x + h, y - h) - f(x - h, y + h) + f(x - h, y - h)) / (4*h*h);
24     return d;
25 }
26 
27 // Gaussian curvature for Monge patch z = f(x,y):
28 // K = (f_xx f_yy - f_xy^2) / (1 + f_x^2 + f_y^2)^2
29 static double gaussian_curvature(double x, double y, double h = 1e-4) {
30     Derivs d = derivatives(x, y, h);
31     double num = d.fxx * d.fyy - d.fxy * d.fxy;
32     double den = pow(1.0 + d.fx*d.fx + d.fy*d.fy, 2);
33     return num / den;
34 }
35 
36 int main() {
37     cout.setf(std::ios::fixed); cout<<setprecision(6);
38 
39     // Evaluate curvature at some points
40     vector<pair<double,double>> pts = { {0,0}, {0.5,0.5}, {1.0,0.0}, {0.0,1.0} };
41     for (auto [x,y] : pts) {
42         double K = gaussian_curvature(x,y);
43         cout << "Point (" << x << ", " << y << ") -> K = " << K << "\n";
44     }
45 
46     // Optional: sample on a grid and compute min/max curvature
47     double xmin=-1, xmax=1, ymin=-1, ymax=1; int nx=41, ny=41;
48     double Kmin=1e100, Kmax=-1e100;
49     for(int i=0;i<nx;i++){
50         double x = xmin + (xmax-xmin)*i/(nx-1);
51         for(int j=0;j<ny;j++){
52             double y = ymin + (ymax-ymin)*j/(ny-1);
53             double K = gaussian_curvature(x,y);
54             Kmin = min(Kmin, K);
55             Kmax = max(Kmax, K);
56         }
57     }
58     cout << "Grid K range: [" << Kmin << ", " << Kmax << "]\n";
59     return 0;
60 }
61

This program computes Gaussian curvature for a height field z = f(x, y) using central finite differences. It estimates first and second derivatives of f and plugs them into the Monge patch formula K = (f_xx f_yy − f_xy^2)/(1 + f_x^2 + f_y^2)^2. The example function is a saddle, so K should be negative or near zero depending on location. The code demonstrates per-point O(1) computation and can be extended to grids.

Time: O(1) per point (O(nm) for an n×m grid)Space: O(1) per point (O(nm) if storing a grid)

2D curvature from a metric: Christoffel, Riemann, Gaussian K, and scalar curvature on S^2

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 // Small 2x2 linear algebra helpers
5 struct Mat2 {
6     double a11, a12, a21, a22; // row-major
7 };
8 
9 static Mat2 mat_inv(const Mat2 &A) {
10     double det = A.a11*A.a22 - A.a12*A.a21;
11     return { A.a22/det, -A.a12/det, -A.a21/det, A.a11/det };
12 }
13 static double mat_det(const Mat2 &A){ return A.a11*A.a22 - A.a12*A.a21; }
14 
15 // Metric on a 2-sphere of radius R in (theta, phi) coordinates:
16 // g = diag(R^2, R^2 sin^2 theta)
17 struct SphereMetric {
18     double R;
19     Mat2 g(double theta, double phi) const {
20         (void)phi; // metric is independent of phi
21         double s = sin(theta);
22         return { R*R, 0.0, 0.0, R*R*s*s };
23     }
24 };
25 
26 // Compute Christoffel symbols Gamma^k_{ij} numerically from metric using central differences
27 // i,j,k in {0,1} correspond to {theta, phi}
28 struct Connection {
29     array<array<array<double,2>,2>,2> Gamma; // Gamma[k][i][j]
30 };
31 
32 static Connection christoffel_at(const SphereMetric &M, double th, double ph, double h=1e-6){
33     // Metric and inverse
34     Mat2 g = M.g(th, ph);
35     Mat2 gi = mat_inv(g);
36 
37     // Partial derivatives of g along theta and phi via central differences
38     auto g_dir = [&](int dir, double d){ // dir=0:theta, 1:phi
39         double th2 = th + (dir==0? d:0.0);
40         double ph2 = ph + (dir==1? d:0.0);
41         return M.g(th2, ph2);
42     };
43     auto dg = [&](int dir){
44         Mat2 gp = g_dir(dir, +h), gm = g_dir(dir, -h);
45         return Mat2{ (gp.a11-gm.a11)/(2*h), (gp.a12-gm.a12)/(2*h), (gp.a21-gm.a21)/(2*h), (gp.a22-gm.a22)/(2*h) };
46     };
47     Mat2 dg_th = dg(0), dg_ph = dg(1);
48 
49     auto get = [&](const Mat2 &A, int r, int c){ return (r==0 && c==0)?A.a11 : (r==0 && c==1)?A.a12 : (r==1 && c==0)?A.a21 : A.a22; };
50 
51     Connection C{};
52     for(int i=0;i<2;i++) for(int j=0;j<2;j++) for(int k=0;k<2;k++){
53         double sum = 0.0;
54         for(int ell=0; ell<2; ++ell){
55             // partial_i g_{ell j}
56             double d_i_g_ellj = (i==0? get(dg_th, ell, j) : get(dg_ph, ell, j));
57             // partial_j g_{ell i}
58             double d_j_g_elli = (j==0? get(dg_th, ell, i) : get(dg_ph, ell, i));
59             // partial_ell g_{ij}
60             // compute derivative along ell by finite difference on that coordinate
61             Mat2 dg_ell = (ell==0? dg_th : dg_ph);
62             double d_ell_g_ij = get(dg_ell, i, j);
63             // g^{k ell}
64             double ginv_kell = (k==0 && ell==0)? gi.a11 : (k==0 && ell==1)? gi.a12 : (k==1 && ell==0)? gi.a21 : gi.a22;
65             sum += 0.5 * ginv_kell * ( d_i_g_ellj + d_j_g_elli - d_ell_g_ij );
66         }
67         C.Gamma[k][i][j] = sum;
68     }
69     return C;
70 }
71 
72 // Numerical derivative of Gamma along a coordinate direction
73 static Connection dGamma_dir(const SphereMetric &M, double th, double ph, int dir, double h=1e-6){
74     Connection Cp = christoffel_at(M, th + (dir==0? h:0.0), ph + (dir==1? h:0.0), h);
75     Connection Cm = christoffel_at(M, th - (dir==0? h:0.0), ph - (dir==1? h:0.0), h);
76     Connection D{};
77     for(int k=0;k<2;k++) for(int i=0;i<2;i++) for(int j=0;j<2;j++){
78         D.Gamma[k][i][j] = (Cp.Gamma[k][i][j] - Cm.Gamma[k][i][j]) / (2*h);
79     }
80     return D;
81 }
82 
83 int main(){
84     cout.setf(std::ios::fixed); cout<<setprecision(8);
85 
86     SphereMetric M{1.0}; // unit sphere, true K = 1
87     double th = 1.0;     // avoid poles (0 or pi) due to coordinate singularity
88     double ph = 0.7;
89     double h = 1e-5;
90 
91     // Metric, inverse, and determinant at point
92     Mat2 g = M.g(th, ph);
93     Mat2 gi = mat_inv(g);
94     double detg = mat_det(g);
95 
96     // Christoffel and their derivatives
97     Connection Gam = christoffel_at(M, th, ph, h);
98     Connection dGam_th = dGamma_dir(M, th, ph, 0, h);
99     Connection dGam_ph = dGamma_dir(M, th, ph, 1, h);
100 
101     // R^l_{ijk} = d_j Gamma^l_{ik} - d_i Gamma^l_{jk} + Gamma^l_{jm} Gamma^m_{ik} - Gamma^l_{im} Gamma^m_{jk}
102     // We'll assemble R^l_{i j k} for i,j,k in {0,1}
103     double R_up[2][2][2][2] = {};
104     auto dGam = [&](int dir, int l, int i, int k){ return (dir==0? dGam_th.Gamma[l][i][k] : dGam_ph.Gamma[l][i][k]); };
105 
106     for(int l=0;l<2;l++) for(int i=0;i<2;i++) for(int j=0;j<2;j++) for(int k=0;k<2;k++){
107         double term = dGam(j, l, i, k) - dGam(i, l, j, k);
108         double quad = 0.0;
109         for(int m=0;m<2;m++){
110             quad += Gam.Gamma[l][j][m]*Gam.Gamma[m][i][k] - Gam.Gamma[l][i][m]*Gam.Gamma[m][j][k];
111         }
112         R_up[l][i][j][k] = term + quad;
113     }
114 
115     // Lower the first index: R_{lijk} = g_{l m} R^m_{ijk}
116     double R_low[2][2][2][2] = {};
117     auto g_rc = [&](int r, int c){ return (r==0&&c==0)? g.a11 : (r==0&&c==1)? g.a12 : (r==1&&c==0)? g.a21 : g.a22; };
118     for(int l=0;l<2;l++) for(int i=0;i<2;i++) for(int j=0;j<2;j++) for(int k=0;k<2;k++){
119         double sum = 0.0;
120         for(int m=0;m<2;m++) sum += g_rc(l,m) * R_up[m][i][j][k];
121         R_low[l][i][j][k] = sum;
122     }
123 
124     // In 2D, Gaussian curvature K = R_{0101} / det(g) (using 0-based indices for 1,2)
125     double R0101 = R_low[0][1][0][1];
126     double K = R0101 / detg;
127     double scalar_R = 2.0 * K; // in 2D
128 
129     cout << "det(g) = " << detg << "\n";
130     cout << "K (numeric) = " << K << " (expected 1.0)\n";
131     cout << "Scalar curvature R = " << scalar_R << " (expected 2.0)\n";
132 
133     return 0;
134 }
135

This program computes curvature on the 2-sphere using only the metric in spherical coordinates. It numerically differentiates the metric to get Christoffel symbols, differentiates them again to build the Riemann tensor, and then extracts Gaussian curvature using K = R_{1212}/det(g). In 2D the scalar curvature equals 2K. Results agree with theory (K = 1 for the unit sphere) away from coordinate singularities.

Time: O(1) per evaluation (constant-size 2×2 tensors and a fixed number of finite differences)Space: O(1)

Discrete Gaussian curvature on a triangle mesh via angle deficit

1 #include <bits/stdc++.h>
2 using namespace std;
3 
4 struct Vec3{ double x,y,z; };
5 static Vec3 operator+(const Vec3&a,const Vec3&b){ return {a.x+b.x,a.y+b.y,a.z+b.z}; }
6 static Vec3 operator-(const Vec3&a,const Vec3&b){ return {a.x-b.x,a.y-b.y,a.z-b.z}; }
7 static Vec3 operator*(const Vec3&a,double s){ return {a.x*s,a.y*s,a.z*s}; }
8 static double dot(const Vec3&a,const Vec3&b){ return a.x*b.x+a.y*b.y+a.z*b.z; }
9 static Vec3 cross(const Vec3&a,const Vec3&b){ return {a.y*b.z-a.z*b.y, a.z*b.x-a.x*b.z, a.x*b.y-a.y*b.x}; }
10 static double norm(const Vec3&a){ return sqrt(dot(a,a)); }
11 
12 struct Face{ int a,b,c; };
13 
14 static double triangle_area(const Vec3&A,const Vec3&B,const Vec3&C){ return 0.5*norm(cross(B-A,C-A)); }
15 static double angle_at(const Vec3&A,const Vec3&B,const Vec3&C){
16     // angle at A in triangle ABC
17     Vec3 u = B - A, v = C - A;
18     double lu = norm(u), lv = norm(v);
19     double c = dot(u,v)/(max(1e-15, lu*lv));
20     c = max(-1.0, min(1.0, c));
21     return acos(c);
22 }
23 
24 int main(){
25     cout.setf(std::ios::fixed); cout<<setprecision(6);
26 
27     // Build an octahedron (unit sphere approximation) with 8 faces, 6 vertices
28     vector<Vec3> V = {
29         { 1,0,0}, {-1,0,0}, {0, 1,0}, {0,-1,0}, {0,0, 1}, {0,0,-1}
30     };
31     vector<Face> F = {
32         {0,4,2},{2,4,1},{1,4,3},{3,4,0}, // top pyramid
33         {2,5,0},{1,5,2},{3,5,1},{0,5,3}  // bottom pyramid
34     };
35 
36     int n = (int)V.size();
37     vector<double> angleSum(n, 0.0), areaSum(n, 0.0);
38 
39     for(const auto &t : F){
40         const Vec3 &A = V[t.a], &B = V[t.b], &C = V[t.c];
41         double area = triangle_area(A,B,C);
42         double aA = angle_at(A,B,C);
43         double aB = angle_at(B,C,A);
44         double aC = angle_at(C,A,B);
45         angleSum[t.a] += aA;
46         angleSum[t.b] += aB;
47         angleSum[t.c] += aC;
48         // Barycentric area: 1/3 of face area to each vertex
49         areaSum[t.a] += area/3.0;
50         areaSum[t.b] += area/3.0;
51         areaSum[t.c] += area/3.0;
52     }
53 
54     const double TWO_PI = 2.0 * M_PI;
55     for(int i=0;i<n;i++){
56         double deficit = TWO_PI - angleSum[i];
57         double K = deficit / max(1e-15, areaSum[i]);
58         cout << "Vertex " << i << ": angleSum=" << angleSum[i]
59              << ", area=" << areaSum[i] << ", K~=" << K << "\n";
60     }
61 
62     return 0;
63 }
64

This program estimates discrete Gaussian curvature at mesh vertices using the angle-deficit formula K ≈ (2π − sum of incident angles)/A_v, with A_v taken as one-third of incident triangle areas. On the octahedron (a coarse sphere), curvature is positive at all vertices, reflecting spherical geometry. Finer meshes and better area schemes (e.g., mixed Voronoi) improve accuracy.

Time: O(F) to accumulate angles and areas (three corners per face)Space: O(V + F) for mesh storage and O(V) for curvature arrays

1	#include <bits/stdc++.h>
2	using namespace std;
3
4	// Example height function: a saddle z = (x^2 - y^2)/2
5	static double f(double x, double y) {
6	return 0.5 * (xx - yy);
7	}
8
9	// Central finite differences for first and second derivatives
10	struct Derivs {
11	double fx, fy, fxx, fyy, fxy;
12	};
13
14	Derivs derivatives(double x, double y, double h = 1e-4) {
15	Derivs d{};
16	// First derivatives
17	d.fx = (f(x + h, y) - f(x - h, y)) / (2*h);
18	d.fy = (f(x, y + h) - f(x, y - h)) / (2*h);
19	// Second derivatives
20	d.fxx = (f(x + h, y) - 2f(x, y) + f(x - h, y)) / (hh);
21	d.fyy = (f(x, y + h) - 2f(x, y) + f(x, y - h)) / (hh);
22	// Mixed derivative via central differences
23	d.fxy = (f(x + h, y + h) - f(x + h, y - h) - f(x - h, y + h) + f(x - h, y - h)) / (4hh);
24	return d;
25	}
26
27	// Gaussian curvature for Monge patch z = f(x,y):
28	// K = (f_xx f_yy - f_xy^2) / (1 + f_x^2 + f_y^2)^2
29	static double gaussian_curvature(double x, double y, double h = 1e-4) {
30	Derivs d = derivatives(x, y, h);
31	double num = d.fxx * d.fyy - d.fxy * d.fxy;
32	double den = pow(1.0 + d.fxd.fx + d.fyd.fy, 2);
33	return num / den;
34	}
35
36	int main() {
37	cout.setf(std::ios::fixed); cout<<setprecision(6);
38
39	// Evaluate curvature at some points
40	vector<pair<double,double>> pts = { {0,0}, {0.5,0.5}, {1.0,0.0}, {0.0,1.0} };
41	for (auto [x,y] : pts) {
42	double K = gaussian_curvature(x,y);
43	cout << "Point (" << x << ", " << y << ") -> K = " << K << "\n";
44	}
45
46	// Optional: sample on a grid and compute min/max curvature
47	double xmin=-1, xmax=1, ymin=-1, ymax=1; int nx=41, ny=41;
48	double Kmin=1e100, Kmax=-1e100;
49	for(int i=0;i<nx;i++){
50	double x = xmin + (xmax-xmin)*i/(nx-1);
51	for(int j=0;j<ny;j++){
52	double y = ymin + (ymax-ymin)*j/(ny-1);
53	double K = gaussian_curvature(x,y);
54	Kmin = min(Kmin, K);
55	Kmax = max(Kmax, K);
56	}
57	}
58	cout << "Grid K range: [" << Kmin << ", " << Kmax << "]\n";
59	return 0;
60	}
61

1	#include <bits/stdc++.h>
2	using namespace std;
3
4	// Small 2x2 linear algebra helpers
5	struct Mat2 {
6	double a11, a12, a21, a22; // row-major
7	};
8
9	static Mat2 mat_inv(const Mat2 &A) {
10	double det = A.a11A.a22 - A.a12A.a21;
11	return { A.a22/det, -A.a12/det, -A.a21/det, A.a11/det };
12	}
13	static double mat_det(const Mat2 &A){ return A.a11A.a22 - A.a12A.a21; }
14
15	// Metric on a 2-sphere of radius R in (theta, phi) coordinates:
16	// g = diag(R^2, R^2 sin^2 theta)
17	struct SphereMetric {
18	double R;
19	Mat2 g(double theta, double phi) const {
20	(void)phi; // metric is independent of phi
21	double s = sin(theta);
22	return { RR, 0.0, 0.0, RRss };
23	}
24	};
25
26	// Compute Christoffel symbols Gamma^k_{ij} numerically from metric using central differences
27	// i,j,k in {0,1} correspond to {theta, phi}
28	struct Connection {
29	array<array<array<double,2>,2>,2> Gamma; // Gamma[k][i][j]
30	};
31
32	static Connection christoffel_at(const SphereMetric &M, double th, double ph, double h=1e-6){
33	// Metric and inverse
34	Mat2 g = M.g(th, ph);
35	Mat2 gi = mat_inv(g);
36
37	// Partial derivatives of g along theta and phi via central differences
38	auto g_dir = [&](int dir, double d){ // dir=0:theta, 1:phi
39	double th2 = th + (dir==0? d:0.0);
40	double ph2 = ph + (dir==1? d:0.0);
41	return M.g(th2, ph2);
42	};
43	auto dg = [&](int dir){
44	Mat2 gp = g_dir(dir, +h), gm = g_dir(dir, -h);
45	return Mat2{ (gp.a11-gm.a11)/(2h), (gp.a12-gm.a12)/(2h), (gp.a21-gm.a21)/(2h), (gp.a22-gm.a22)/(2h) };
46	};
47	Mat2 dg_th = dg(0), dg_ph = dg(1);
48
49	auto get = [&](const Mat2 &A, int r, int c){ return (r==0 && c==0)?A.a11 : (r==0 && c==1)?A.a12 : (r==1 && c==0)?A.a21 : A.a22; };
50
51	Connection C{};
52	for(int i=0;i<2;i++) for(int j=0;j<2;j++) for(int k=0;k<2;k++){
53	double sum = 0.0;
54	for(int ell=0; ell<2; ++ell){
55	// partial_i g_{ell j}
56	double d_i_g_ellj = (i==0? get(dg_th, ell, j) : get(dg_ph, ell, j));
57	// partial_j g_{ell i}
58	double d_j_g_elli = (j==0? get(dg_th, ell, i) : get(dg_ph, ell, i));
59	// partial_ell g_{ij}
60	// compute derivative along ell by finite difference on that coordinate
61	Mat2 dg_ell = (ell==0? dg_th : dg_ph);
62	double d_ell_g_ij = get(dg_ell, i, j);
63	// g^{k ell}
64	double ginv_kell = (k==0 && ell==0)? gi.a11 : (k==0 && ell==1)? gi.a12 : (k==1 && ell==0)? gi.a21 : gi.a22;
65	sum += 0.5 * ginv_kell * ( d_i_g_ellj + d_j_g_elli - d_ell_g_ij );
66	}
67	C.Gamma[k][i][j] = sum;
68	}
69	return C;
70	}
71
72	// Numerical derivative of Gamma along a coordinate direction
73	static Connection dGamma_dir(const SphereMetric &M, double th, double ph, int dir, double h=1e-6){
74	Connection Cp = christoffel_at(M, th + (dir==0? h:0.0), ph + (dir==1? h:0.0), h);
75	Connection Cm = christoffel_at(M, th - (dir==0? h:0.0), ph - (dir==1? h:0.0), h);
76	Connection D{};
77	for(int k=0;k<2;k++) for(int i=0;i<2;i++) for(int j=0;j<2;j++){
78	D.Gamma[k][i][j] = (Cp.Gamma[k][i][j] - Cm.Gamma[k][i][j]) / (2*h);
79	}
80	return D;
81	}
82
83	int main(){
84	cout.setf(std::ios::fixed); cout<<setprecision(8);
85
86	SphereMetric M{1.0}; // unit sphere, true K = 1
87	double th = 1.0; // avoid poles (0 or pi) due to coordinate singularity
88	double ph = 0.7;
89	double h = 1e-5;
90
91	// Metric, inverse, and determinant at point
92	Mat2 g = M.g(th, ph);
93	Mat2 gi = mat_inv(g);
94	double detg = mat_det(g);
95
96	// Christoffel and their derivatives
97	Connection Gam = christoffel_at(M, th, ph, h);
98	Connection dGam_th = dGamma_dir(M, th, ph, 0, h);
99	Connection dGam_ph = dGamma_dir(M, th, ph, 1, h);
100
101	// R^l_{ijk} = d_j Gamma^l_{ik} - d_i Gamma^l_{jk} + Gamma^l_{jm} Gamma^m_{ik} - Gamma^l_{im} Gamma^m_{jk}
102	// We'll assemble R^l_{i j k} for i,j,k in {0,1}
103	double R_up[2][2][2][2] = {};
104	auto dGam = [&](int dir, int l, int i, int k){ return (dir==0? dGam_th.Gamma[l][i][k] : dGam_ph.Gamma[l][i][k]); };
105
106	for(int l=0;l<2;l++) for(int i=0;i<2;i++) for(int j=0;j<2;j++) for(int k=0;k<2;k++){
107	double term = dGam(j, l, i, k) - dGam(i, l, j, k);
108	double quad = 0.0;
109	for(int m=0;m<2;m++){
110	quad += Gam.Gamma[l][j][m]Gam.Gamma[m][i][k] - Gam.Gamma[l][i][m]Gam.Gamma[m][j][k];
111	}
112	R_up[l][i][j][k] = term + quad;
113	}
114
115	// Lower the first index: R_{lijk} = g_{l m} R^m_{ijk}
116	double R_low[2][2][2][2] = {};
117	auto g_rc = [&](int r, int c){ return (r==0&&c==0)? g.a11 : (r==0&&c==1)? g.a12 : (r==1&&c==0)? g.a21 : g.a22; };
118	for(int l=0;l<2;l++) for(int i=0;i<2;i++) for(int j=0;j<2;j++) for(int k=0;k<2;k++){
119	double sum = 0.0;
120	for(int m=0;m<2;m++) sum += g_rc(l,m) * R_up[m][i][j][k];
121	R_low[l][i][j][k] = sum;
122	}
123
124	// In 2D, Gaussian curvature K = R_{0101} / det(g) (using 0-based indices for 1,2)
125	double R0101 = R_low[0][1][0][1];
126	double K = R0101 / detg;
127	double scalar_R = 2.0 * K; // in 2D
128
129	cout << "det(g) = " << detg << "\n";
130	cout << "K (numeric) = " << K << " (expected 1.0)\n";
131	cout << "Scalar curvature R = " << scalar_R << " (expected 2.0)\n";
132
133	return 0;
134	}
135

1	#include <bits/stdc++.h>
2	using namespace std;
3
4	struct Vec3{ double x,y,z; };
5	static Vec3 operator+(const Vec3&a,const Vec3&b){ return {a.x+b.x,a.y+b.y,a.z+b.z}; }
6	static Vec3 operator-(const Vec3&a,const Vec3&b){ return {a.x-b.x,a.y-b.y,a.z-b.z}; }
7	static Vec3 operator(const Vec3&a,double s){ return {a.xs,a.ys,a.zs}; }
8	static double dot(const Vec3&a,const Vec3&b){ return a.xb.x+a.yb.y+a.z*b.z; }
9	static Vec3 cross(const Vec3&a,const Vec3&b){ return {a.yb.z-a.zb.y, a.zb.x-a.xb.z, a.xb.y-a.yb.x}; }
10	static double norm(const Vec3&a){ return sqrt(dot(a,a)); }
11
12	struct Face{ int a,b,c; };
13
14	static double triangle_area(const Vec3&A,const Vec3&B,const Vec3&C){ return 0.5*norm(cross(B-A,C-A)); }
15	static double angle_at(const Vec3&A,const Vec3&B,const Vec3&C){
16	// angle at A in triangle ABC
17	Vec3 u = B - A, v = C - A;
18	double lu = norm(u), lv = norm(v);
19	double c = dot(u,v)/(max(1e-15, lu*lv));
20	c = max(-1.0, min(1.0, c));
21	return acos(c);
22	}
23
24	int main(){
25	cout.setf(std::ios::fixed); cout<<setprecision(6);
26
27	// Build an octahedron (unit sphere approximation) with 8 faces, 6 vertices
28	vector<Vec3> V = {
29	{ 1,0,0}, {-1,0,0}, {0, 1,0}, {0,-1,0}, {0,0, 1}, {0,0,-1}
30	};
31	vector<Face> F = {
32	{0,4,2},{2,4,1},{1,4,3},{3,4,0}, // top pyramid
33	{2,5,0},{1,5,2},{3,5,1},{0,5,3} // bottom pyramid
34	};
35
36	int n = (int)V.size();
37	vector<double> angleSum(n, 0.0), areaSum(n, 0.0);
38
39	for(const auto &t : F){
40	const Vec3 &A = V[t.a], &B = V[t.b], &C = V[t.c];
41	double area = triangle_area(A,B,C);
42	double aA = angle_at(A,B,C);
43	double aB = angle_at(B,C,A);
44	double aC = angle_at(C,A,B);
45	angleSum[t.a] += aA;
46	angleSum[t.b] += aB;
47	angleSum[t.c] += aC;
48	// Barycentric area: 1/3 of face area to each vertex
49	areaSum[t.a] += area/3.0;
50	areaSum[t.b] += area/3.0;
51	areaSum[t.c] += area/3.0;
52	}
53
54	const double TWO_PI = 2.0 * M_PI;
55	for(int i=0;i<n;i++){
56	double deficit = TWO_PI - angleSum[i];
57	double K = deficit / max(1e-15, areaSum[i]);
58	cout << "Vertex " << i << ": angleSum=" << angleSum[i]
59	<< ", area=" << areaSum[i] << ", K~=" << K << "\n";
60	}
61
62	return 0;
63	}
64