feat(mesh): mesh quality optimisation (smooth, aspect, edge flips) (P33)

[MAfarrag]

Context

Phase 3's P26 (#40) exports meshes to Gmsh, Iber, HEC-RAS, etc. Raw exports from a structured DEM produce
meshes with very poor quality near breaklines and building footprints — sliver triangles, near-zero areas,
and inverted elements. Hydrodynamic solvers either reject these meshes outright or run them with severe
time-step penalties. Phase 4 adds a post-export quality-optimisation pass.

Problem / Current Behaviour

After P26 ships, mesh files load in Gmsh but visibly contain slivers, obtuse triangles, and non-conforming
boundary nodes. Production workflows require manually opening Gmsh and running its optimiser before
hydraulic simulation — friction that should be in-process.

Affected locations

File	Symbol	Notes
src/digitalrivers/mesh.py	Mesh.laplacian_smooth(...)	shipped — 1-ring centroid relaxation.
src/digitalrivers/mesh.py	Mesh.aspect_ratios()	shipped — per-triangle quality metric.
src/digitalrivers/mesh.py	(edge swap / Delaunay refinement)	pending.
src/digitalrivers/mesh.py	(node insertion at low-quality triangles)	pending.
src/digitalrivers/mesh.py	(boundary projection back onto breaklines)	pending.
digitalrivers.export.mesh.optimise_mesh(...)	(top-level wrapper)	pending.

Motivation Example

from digitalrivers import DEM
from digitalrivers.mesh import Mesh

dem  = DEM.read_file("conditioned.tif")
mesh = dem.to_mesh()                                # exists via P26 (#40)

# Shipped — node smoothing + quality reporting.
smoothed = mesh.laplacian_smooth(max_iter=50, tolerance=1e-4)
print(smoothed.aspect_ratios().max())

# Pending — full optimisation pipeline.
# from digitalrivers.export.mesh import optimise_mesh
# optimised, quality = optimise_mesh(
#     smoothed,
#     target_quality=0.3,
#     passes=["swap", "laplacian"],
# )
# print(quality.converged, quality.skewness_max)

Proposed Solution

from digitalrivers.export.mesh import optimise_mesh, MeshQuality

def optimise_mesh(
    mesh: "meshio.Mesh",
    *,
    target_quality: float = 0.3,            # max skewness allowed
    max_iter: int = 100,
    passes: list["Literal['laplacian','swap','insert','project']"] = ["swap", "laplacian"],
    boundaries: "meshio.Mesh | None" = None,
) -> tuple["meshio.Mesh", "MeshQuality"]: ...


@dataclass
class MeshQuality:
    n_elements: int
    n_nodes: int
    aspect_ratio_max: float
    skewness_max: float
    area_ratio: float
    edge_length_ratio: float
    converged: bool
    iterations: int

Four optimisation passes:

• Laplacian smoothing: move each interior node to centroid of 1-ring neighbours; iterate until movement
  below tolerance. Preserves boundaries.
• Edge swapping (Delaunay refinement): for each interior edge shared by two triangles, swap if it reduces
  the maximum angle. Lawson convergence proof for local-Delaunay.
• Node insertion: subdivide triangles whose quality is below threshold; new node at the circumcentre
  (Bowyer-Watson incremental Delaunay).
• Boundary projection: ensure boundary nodes stay on the input breakline / DEM-edge polylines.

Quality metric: Liu-Joe form of the radius ratio, q = 4·√3·A / (l1² + l2² + l3²) — algebraically equivalent
to 2·r_in / r_out but cheaper to compute.

References

• George P. L., Borouchaki H. (1998). Delaunay Triangulation and Meshing. Hermès.
• Freitag L. A., Plassmann P. (2000). "Local optimization-based simplicial mesh untangling and improvement."
  Int. J. Num. Meth. Eng. 49.
• Shewchuk J. R. (2002). "What is a good linear finite element?" — radius-ratio quality metric.
• Gmsh gmsh.model.mesh.optimize(method, ...) and pygmsh._optimize.
• meshpy for Triangle/TetGen integration.

Out of Scope

• 3D / volumetric mesh optimisation (tetrahedral) — Phase 5+.
• Higher-order curved elements — out-of-package.

Effort Estimate

Size: L (5 days for core four passes + 2 days for boundary-projection edge cases)
Rationale: Laplacian smoothing + quality metrics are 1 day. Edge swap + node insertion are 2–3 days
including correctness tests. Boundary projection edge cases are the slow part.

Definition of Done

• Mesh.laplacian_smooth(...) ships and is exercised by tests.
• Mesh.aspect_ratios() ships and reports per-triangle quality.
• Pending: optimise_mesh(...) wrapper with the four-pass pipeline.
• Pending: Synthetic mesh with a sliver triangle: after optimise_mesh(passes=["swap"]), the sliver
  is eliminated.
• Pending: Laplacian smoothing preserves boundary nodes (fixed-node test).
• Pending: Raw mesh from a 1 m DEM with breaklines: target_quality=0.3 achievable in <= 100
  iterations.
• Pending: Output mesh loads cleanly in Gmsh GUI (eye test, pytest.mark.visual).
• Pending: No degenerate elements after optimisation (mesh.cells[...].area > 1e-9).

Status: Partial on branch feat/phase-4. Shipped: laplacian_smooth + aspect_ratios. Pending: edge
swap, node insertion, boundary projection, optimise_mesh wrapper. See tests/test_mesh.py.