# Discrete Calculus on Simplicial Meshes

## Overview

This module implements differential operators (gradient, divergence, curl,
Laplacian) for simplicial meshes using two complementary approaches:

1. **Discrete Exterior Calculus (DEC)** - Rigorous differential geometry
   framework based on Desbrun et al. (2005) and Hirani (2003)
2. **Weighted Least-Squares (LSQ)** - Practical CFD/FEM approach for general
   use cases

---

## Discrete Exterior Calculus (DEC)

DEC provides a mathematically rigorous framework where discrete operators
satisfy exact discrete versions of continuous theorems (Stokes, Gauss-Bonnet,
etc.).

### Core DEC Operators

#### Laplace-Beltrami Operator

```python
from physicsnemo.mesh.calculus.laplacian import compute_laplacian_points_dec

# Intrinsic Laplacian: Δf(v) = -(1/|⋆v|) Σ (|⋆e|/|e|)(f_neighbor - f_v)
laplacian = compute_laplacian_points_dec(mesh, scalar_field)
```

**Properties**:

- Uses cotangent weights derived from the FEM stiffness matrix (see below)
- In 2D: equivalent to `(1/2)(cot α + cot β)` (Meyer Eq. 5)
- In 1D: gives `1/|edge|` (standard finite-difference second derivative)
- In 3D+: gives exact dihedral-angle-based weights via Gram matrix inverse
- Normalized by circumcentric dual volumes (Voronoi cells)
- Exact for linear functions at interior vertices
- Works on manifolds of any dimension embedded in any ambient space

**Cotangent weights via FEM stiffness matrix** (n-dimensional):

For an n-simplex with vertices v₀, ..., vₙ, the cotangent weight for
edge (i, j) is:

```text
w_ij = -|σ| × (∇λ_i · ∇λ_j)
```

where λ_i are barycentric coordinate functions and |σ| is the cell volume.
The gradient dot products are computed from the inverse Gram matrix:

```text
E = [v₁-v₀, ..., vₙ-v₀]    (edge matrix, n × d)
G = E @ E^T                   (Gram matrix, n × n)
∇λ_k · ∇λ_l = (G⁻¹)_{k-1,l-1}   for k, l ≥ 1
```

This generalizes the classical 2D cotangent formula to arbitrary dimensions.

**Reference**: Hirani (2003) Eq. 6.4.2, Meyer et al. (2003) Eq. 8

#### Exterior Derivative

```python
from physicsnemo.mesh.calculus._exterior_derivative import (
    exterior_derivative_0,
    exterior_derivative_1,
)

# d: Ω⁰ → Ω¹ (0-forms to 1-forms)
# df([vi,vj]) = f(vj) - f(vi)
edge_1form, edges = exterior_derivative_0(mesh, vertex_values)

# d: Ω¹ → Ω² (1-forms to 2-forms)
# Circulation around faces
face_2form, faces = exterior_derivative_1(mesh, edge_1form, edges)
```

**Properties**:

- `d ∘ d = 0` (exact by construction)
- Discrete Stokes theorem: `⟨dα, c⟩ = ⟨α, ∂c⟩` (true by definition)

**Reference**: Desbrun et al. (2005) Section 3, Hirani (2003) Chapter 3

#### Hodge Star

```python
from physicsnemo.mesh.calculus._hodge_star import hodge_star_0, hodge_star_1

# ⋆: Ω⁰ → Ωⁿ (vertex values to dual n-cells)
star_f = hodge_star_0(mesh, f)  # ⋆f(⋆v) = f(v) × |⋆v|
```

**Properties**:

- Preserves averages: `⟨α, σ⟩/|σ| = ⟨⋆α, ⋆σ⟩/|⋆σ|`
- `⋆⋆α = (-1)^(k(n-k)) α`
- Uses circumcentric (Voronoi) dual cells, NOT barycentric

**Reference**: Hirani (2003) Def. 4.1.1, Desbrun et al. (2005) Section 4

#### Sharp and Flat Operators

```python
from physicsnemo.mesh.calculus._sharp_flat import sharp, flat

# ♯: Ω¹ → 𝔛 (1-forms to vector fields)
grad_vector = sharp(mesh, df, edges)

# ♭: 𝔛 → Ω¹ (vector fields to 1-forms)
one_form = flat(mesh, vector_field, edges)
```

**Implementation**:

- **Sharp (♯)**: Hirani Eq. 5.8.1 with support volume intersections and
  barycentric gradients
- **Flat (♭)**: PDP-flat (Hirani Section 5.6) using averaged endpoint vectors

**Note**: Sharp and flat are NOT exact inverses in discrete DEC (Hirani
Prop. 5.5.3). This is a fundamental property of the discrete theory, not a bug.

**Reference**: Hirani (2003) Chapter 5

### Gradient via DEC

```python
from physicsnemo.mesh.calculus.gradient import compute_gradient_points_dec

# Computes: grad(f) = ♯(df)
grad_f = compute_gradient_points_dec(mesh, scalar_field)
```

Combines exterior derivative and sharp operator to produce gradient vector
field.

---

## Weighted Least-Squares (LSQ) Methods

LSQ methods provide general-purpose operators that work robustly on arbitrary
meshes.

### Gradient

```python
from physicsnemo.mesh.calculus.gradient import (
    compute_gradient_points_lsq,
    compute_gradient_cells_lsq,
)

# At vertices
grad = compute_gradient_points_lsq(
    mesh,
    scalar_field,
    weight_power=2.0,  # Inverse distance weighting
    intrinsic=False    # Set True for tangent-space gradients on manifolds
)

# At cell centers
grad_cells = compute_gradient_cells_lsq(mesh, cell_values)
```

**Properties**:

- Exact for constant and linear fields
- First-order accurate O(h) for smooth fields
- Supports intrinsic (tangent-space) computation for embedded manifolds
- Works for both scalar and tensor fields

### Divergence

```python
from physicsnemo.mesh.calculus.divergence import compute_divergence_points_lsq

div_v = compute_divergence_points_lsq(mesh, vector_field)
```

Computes `div(v) = ∂vₓ/∂x + ∂vᵧ/∂y + ∂vᵧ/∂z` via component gradients.

### Curl (3D Only)

```python
from physicsnemo.mesh.calculus.curl import compute_curl_points_lsq

curl_v = compute_curl_points_lsq(mesh, vector_field)  # Requires n_spatial_dims = 3
```

Computes curl from antisymmetric part of Jacobian matrix.

---

## Circumcentric Dual Volumes (Voronoi Cells)

### Implementation

```python
from physicsnemo.mesh.geometry.dual_meshes import compute_dual_volumes_0

dual_vols = compute_dual_volumes_0(mesh)  # |⋆v| for each vertex
```

**Algorithm** (dimension-specific):

**1D manifolds (edges)**:

- Each vertex gets half the length of each incident edge
- Exact for piecewise linear 1-manifolds

**2D manifolds (triangles)**:

- **Acute triangles**: Circumcentric Voronoi formula (Meyer Eq. 7)

```text
|⋆v| = (1/8) Σ (||e||² cot(opposite_angle))
```

- **Obtuse triangles**: Mixed area subdivision (Meyer Fig. 4)

```text
If obtuse at vertex: |⋆v| = area(T)/2
Otherwise: |⋆v| = area(T)/4
```

**3D+ manifolds (tetrahedra, etc.)**:

- Barycentric approximation: `|⋆v| = Σ |cell|/(n+1)`
- Note: Rigorous circumcentric dual requires "well-centered" meshes
  (Desbrun 2005)

**Property**: Perfect tiling: `Σ_vertices |⋆v| = |mesh|` (conservation holds
exactly)

**References**:

- Meyer et al. (2003) Sections 3.2-3.4
- Desbrun et al. (2005) lines 286-395
- Hirani (2003) Def. 2.4.5

---

### Known Behavior (Not Bugs)

**div(grad(f)) ≈ Δf but not exactly**:

- In discrete DEC, sharp (♯) and flat (♭) are NOT exact inverses (Hirani
  Prop. 5.5.3)
- Therefore `div(grad(f))` and `Δf` may differ by ~2-3x on coarse meshes
- Both are O(h) accurate, difference → 0 as mesh refines
- This is a fundamental property of discrete exterior calculus

**3D dual volumes use barycentric approximation**:

- Rigorous circumcentric requires "well-centered" meshes (Desbrun 2005)
- Mixed volume formula for obtuse tetrahedra doesn't exist in literature
- Current barycentric approximation is standard practice and works well

---

## API Reference

### High-Level Interface

```python
# Unified interface for derivatives
mesh_with_grad = mesh.compute_point_derivatives(
    keys=['pressure', 'temperature'],
    method='lsq',  # or 'dec' for Laplacian only
    gradient_type='extrinsic',  # or 'intrinsic' for manifolds
    weight_power=2.0,
)

# Access results
grad_p = mesh_with_grad.point_data['pressure_gradient']  # (n_points, n_spatial_dims)
```

### Direct Operator Calls

```python
from physicsnemo.mesh.calculus import (
    compute_gradient_points_lsq,
    compute_divergence_points_lsq,
    compute_curl_points_lsq,
    compute_laplacian_points_dec,
)

# Gradient (LSQ or DEC)
grad = compute_gradient_points_lsq(mesh, f, weight_power=2.0, intrinsic=False)
grad = compute_gradient_points_dec(mesh, f)  # DEC method

# Divergence
div = compute_divergence_points_lsq(mesh, vector_field)

# Curl (3D only)
curl = compute_curl_points_lsq(mesh, vector_field)

# Laplacian (DEC method)
laplacian = compute_laplacian_points_dec(mesh, scalar_field)
```

---

## Performance

All operations are **fully vectorized** (no Python loops over mesh elements):

- **Gradient/Divergence/Curl**: O(n_points × avg_degree)
- **Laplacian**: O(n_edges), very efficient
- **Dual volumes**: O(n_cells), one-time computation with caching

**Memory**: Minimal overhead, intermediate results cached in `TensorDict`

**Scaling**: Designed for massive meshes (100M+ points on GB200-class GPUs)

---

## Module Structure

```text
src/physicsnemo.mesh/calculus/
├── __init__.py                    # Public API
├── derivatives.py                 # High-level interface (compute_point_derivatives)
├── gradient.py                    # Gradient (LSQ + DEC)
├── divergence.py                  # Divergence (LSQ + DEC)
├── curl.py                        # Curl (LSQ, 3D only)
├── laplacian.py                   # Laplace-Beltrami (DEC)
│
├── _exterior_derivative.py        # DEC: exterior derivative d
├── _hodge_star.py                 # DEC: Hodge star ⋆
├── _sharp_flat.py                 # DEC: sharp ♯ and flat ♭
│
├── _lsq_reconstruction.py         # LSQ: gradient reconstruction (ambient space)
└── _lsq_intrinsic.py             # LSQ: intrinsic gradients (tangent space)
```

```text
src/physicsnemo.mesh/geometry/
├── dual_meshes.py                 # Unified dual 0-cell volumes (Voronoi cells)
├── support_volumes.py             # Support volume intersections for DEC
└── interpolation.py               # Barycentric function gradients
```

---

## Usage Examples

### Example 1: Laplace-Beltrami on Curved Surface

```python
import torch
from physicsnemo.mesh.mesh import Mesh
from physicsnemo.mesh.calculus.laplacian import compute_laplacian_points_dec

# Create surface mesh (e.g., sphere, imported mesh, etc.)
mesh = ...  # 2D surface in 3D

# Add scalar field (e.g., temperature distribution)
temperature = mesh.point_data['temperature']

# Compute intrinsic Laplacian
laplacian = compute_laplacian_points_dec(mesh, temperature)

# Use for diffusion: ∂T/∂t = κ Δ T
mesh.point_data['laplacian_T'] = laplacian
```

### Example 2: Gradient on Manifold (Intrinsic)

```python
from physicsnemo.mesh.calculus.gradient import compute_gradient_points_lsq

# Compute gradient in tangent space (for surface in 3D)
grad_intrinsic = compute_gradient_points_lsq(
    mesh,
    scalar_field,
    intrinsic=True,  # Solves in tangent space
)

# Result is guaranteed perpendicular to surface normal
assert torch.allclose(
    (grad_intrinsic * mesh.point_normals).sum(dim=-1),
    torch.zeros(mesh.n_points),
    atol=1e-6
)
```

### Example 3: Vector Calculus Identities

```python
from physicsnemo.mesh.calculus import (
    compute_gradient_points_lsq,
    compute_divergence_points_lsq,
    compute_curl_points_lsq,
)

# Verify curl(grad(f)) = 0
grad_f = compute_gradient_points_lsq(mesh, scalar_field)
curl_grad_f = compute_curl_points_lsq(mesh, grad_f)
assert torch.allclose(curl_grad_f, torch.zeros_like(curl_grad_f), atol=1e-5)

# Verify div(curl(v)) = 0
curl_v = compute_curl_points_lsq(mesh, vector_field)
div_curl_v = compute_divergence_points_lsq(mesh, curl_v)
assert torch.allclose(div_curl_v, torch.zeros_like(div_curl_v), atol=1e-5)
```

---

## Dimension Support

| Operator | 1D | 2D | 3D | nD |
|----------|----|----|----|----|
| Gradient (LSQ) | ✓ | ✓ | ✓ | ✓ |
| Gradient (DEC) | ✓ | ✓ | ✓ | ✓ |
| Divergence | ✓ | ✓ | ✓ | ✓ |
| Curl (LSQ) | - | - | ✓ | - |
| Laplacian (DEC) | ✓ | ✓ | ✓† | - |
| Hodge star | ✓ | ✓ | ✓* | ✓* |

*Uses barycentric approximation for n ≥ 3

†3D Laplacian uses an inverse-edge-length approximation rather than true
dihedral-angle cotangent weights. Accuracy degrades on poorly-shaped tetrahedra.
Not implemented for n > 3.

---

## Choosing Between DEC and LSQ

**Use DEC when**:

- Need mathematically rigorous operators
- Working with differential geometry (curvatures, etc.)
- Require exact discrete theorems (Stokes, Gauss-Bonnet)
- Computing Laplacian on manifolds

**Use LSQ when**:

- Need general-purpose gradient/divergence/curl
- Working with irregular/poor-quality meshes
- Need robust performance on all mesh types
- Computing derivatives of tensor fields

**Both methods**:

- Are first-order accurate O(h)
- Work on irregular meshes
- Are fully vectorized
- Support GPU acceleration

---

## Limitations and Future Work

### Current Limitations

1. **3D+ Dual Volumes**: Uses barycentric approximation (standard practice)
   - Rigorous circumcentric requires "well-centered" meshes
   - Mixed volume for obtuse tets is an open research problem
   - The Laplacian cotangent *weights* are exact for all dimensions (via FEM
     stiffness matrix); only the dual volume *normalization* uses approximation

2. **Sharp/Flat Not Exact Inverses**: `♯ ∘ ♭ ≠ identity` in discrete DEC
   - This is fundamental to discrete theory (Hirani Prop. 5.5.3)
   - Causes `div(grad) ≈ Δ` (not exact)

3. **Boundary Effects**: Cotangent Laplacian assumes complete 1-ring
   neighborhoods
   - Boundary vertices may show artifacts
   - Set `include_boundary=False` in curvature computations

### Future Enhancements

1. **Well-centered mesh detection** for rigorous 3D dual volumes
2. **Additional DEC operators**: wedge product, interior product, Lie
   derivative
3. **Higher-order LSQ** with extended stencils
4. **Convergence analysis**: Verify O(h²) error as mesh refines
5. **Alternative sharp/flat combinations** (DPP-flat, etc.)

---

## Mathematical Foundations

### Discrete Exterior Calculus

- Exterior forms as cochains (Hirani Chapter 3)
- Circumcentric dual complexes (Desbrun Section 2, Hirani Section 2.4)
- Hodge star via volume ratios (Hirani Def. 4.1.1)
- Sharp/flat with support volumes (Hirani Chapter 5)

### Discrete Differential Geometry

- Meyer mixed Voronoi areas for curvature (Meyer Sections 3.2-3.4)
- Cotangent Laplacian for mean curvature (Meyer Eq. 8)
- Angle defect for Gaussian curvature (Meyer Eq. 9)

### Key Theorems Preserved

- Discrete Stokes theorem (exact)
- Gauss-Bonnet theorem (< 0.001% error numerically)
- Conservation of dual volumes (exact)
- Vector calculus identities: `curl ∘ grad = 0`, `div ∘ curl = 0` (exact)

---

## References

1. **Meyer, M., Desbrun, M., Schröder, P., & Barr, A. H.** (2003). "Discrete
   Differential-Geometry Operators for Triangulated 2-Manifolds". *VisMath*.
   - Sections 3.2-3.4: Mixed Voronoi areas
   - Eq. 5: Cotangent weights
   - Eq. 7: Circumcentric Voronoi formula
   - Eq. 8-9: Mean and Gaussian curvature

2. **Desbrun, M., Hirani, A. N., Leok, M., & Marsden, J. E.** (2005).
   "Discrete Exterior Calculus". *arXiv:math/0508341v2*.
   - Section 2: Circumcentric dual complexes
   - Section 3-4: Exterior derivative and Hodge star
   - Lines 268-275: Cotangent weight derivation

3. **Hirani, A. N.** (2003). "Discrete Exterior Calculus". PhD thesis,
   California Institute of Technology.
   - Chapter 5: Sharp and flat operators
   - Eq. 5.8.1: PP-sharp formula
   - Eq. 6.4.2: Laplace-Beltrami
   - Prop. 5.5.1: Support volume intersections
