Add exact-boundary quality meshing: tessellation.fill_interior

[peterdsharpe]

PhysicsNeMo Pull Request

Description

Part 1 of 3 (part 2: #1801, implicit-domain mesher; part 3: #1802, tutorial). No code dependency on #1801 — land in either order.

Adds exact-boundary quality mesh generation: tessellation.fill_interior takes a closed codimension-one boundary Mesh and fills the enclosed interior with quality simplices.

• Mesh-native: a 2D edge Mesh in — loops in any order/orientation; holes, multiple components, and islands-inside-holes resolved automatically by containment — and a volume Mesh out, on the input's device/dtype.
• Exact-boundary contract: input vertices appear bit-identically (leading rows, input order); boundary facets are only ever subdivided, never moved; deterministic (identical inputs give bitwise-identical outputs).
• Guaranteed quality: every output triangle provably meets the requested minimum angle (Ruppert refinement), with an optional cell-size bound and bound-preserving ODT smoothing.
• Provenance (opt-in): with provenance=True, point_data["boundary_marker"] and point_data["source_point"] map output vertices back to input vertices, enabling per-boundary BCs and data propagation. Off by default — the function claims no keys in the user-owned point_data namespace unless asked.
• Engine is a from-scratch, dependency-free constrained Delaunay pipeline (Bowyer–Watson + Sloan segment recovery + even-odd hole removal); citations in the module docstring. The contract is dimension-generic; n=3 raises NotImplementedError pending exact 3D boundary recovery.
• 93 tests (engine invariants, the Mesh-native contract, and regression tests from two external/internal review rounds) and API docs. Tutorial coverage lives in the dedicated mesh-generation tutorial (Add Tutorial 9: mesh generation (stacked on #1800, #1801) #1802) so tutorial 8 stays scoped to I/O. ~90k guaranteed-quality triangles in ~6 s single-core.

Checklist

• I am familiar with the Contributing Guidelines.
• New or existing tests cover these changes.
• The documentation is up to date with these changes.
• The CHANGELOG.md is up to date with these changes.
• An issue is linked to this pull request.
• If I am implementing a new model or modifying any existing model, I have followed the Models Implementation Coding Standards.

Dependencies

No new dependencies.

Review Process

All PRs are reviewed by the PhysicsNeMo team before merging.

Depending on which files are changed, GitHub may automatically assign a maintainer for review.

We are also testing AI-based code review tools (e.g., Greptile), which may add automated comments with a confidence score.
This score reflects the AI’s assessment of merge readiness and is not a qualitative judgment of your work, nor is
it an indication that the PR will be accepted / rejected.

AI-generated feedback should be reviewed critically for usefulness.
You are not required to respond to every AI comment, but they are intended to help both authors and reviewers.
Please react to Greptile comments with 👍 or 👎 to provide feedback on their accuracy.

[copy-pr-bot]

This pull request requires additional validation before any workflows can run on NVIDIA's runners.

Pull request vetters can view their responsibilities here.

Contributors can view more details about this message here.

[greptile-apps]

Greptile Summary

This PR adds tessellation.fill_interior, a from-scratch constrained Delaunay + Ruppert refinement engine that turns a closed codimension-one boundary Mesh into a quality interior volume mesh with bit-identical input vertex preservation and deterministic, provenance-tagged output.

• delaunay.py (1 618 lines) implements the full stack — Bowyer-Watson insertion, Sloan segment recovery, even-odd hole removal, Ruppert refinement, and quality-gated ODT smoothing — in pure Python + NumPy with careful float64 normalization and an explicit robustness model.
• fill_interior.py wraps the engine with a Mesh-native API: loop extraction, O(n²) containment-depth grouping (any nesting supported), multi-component meshing, and a reordering pass that puts all inherited input vertices first with source_point provenance.
• Tests (58 total across two files) cover algorithmic invariants, adversarial geometry, bitwise determinism, dtype/device round-trips, and boundary-case validation; ASV benchmarks are included.

Important Files Changed

Filename	Overview
physicsnemo/mesh/tessellation/delaunay.py	New 1618-line from-scratch constrained Delaunay + Ruppert refinement engine; algorithm is well-structured and all major paths (Bowyer-Watson, Sloan segment recovery, even-odd hole removal, ODT smoothing) are correctly implemented.
physicsnemo/mesh/tessellation/fill_interior.py	Mesh-native entry point for fill_interior; loop extraction, containment grouping, and source_point provenance reordering are correct, but missing guard for the empty-cells case and first-vertex containment sampling in _group_components could confuse edge cases.
test/mesh/tessellation/test_delaunay.py	Comprehensive 782-line test suite covering algorithmic guarantees, adversarial geometry, determinism, and structural validity; scipy is imported lazily inside one test function rather than at the module level.
test/mesh/tessellation/test_fill_interior.py	208-line Mesh-native contract tests covering provenance, nesting, dtype/device round-trip, orientation-free handling, and validation errors; well matched to the public API contract.
benchmarks/physicsnemo/mesh/tessellation/benchmark_delaunay.py	ASV benchmark suite for delaunay_mesh_2d; straightforward and follows benchmark conventions.
physicsnemo/mesh/tessellation/init.py	Correctly exports fill_interior and polygon_interior_point alongside the existing triangulate; module docstring updated.

_{Reviews (1): Last reviewed commit: "[pre-commit.ci] auto fixes from pre-comm..." | Re-trigger Greptile}

[melo-gonzo]

Hey @peterdsharpe, this is super cool! I worked through the jupyter notebook, which was very helpful. You could consider purging out the metadata for execution counts, times, and some cached paths.

With the help of some agentic investigation, there are three things that seem worth flagging for this initial PR. Since I'm not an expert in this code and meshing, I'll lay out some cases here that were flagged by the investigation. It would be great to get a read on if these issues would be meaningful in real-world usage. 2 and 3 seem valid, while 1 seems like a more general issue with 'reflex corners' and meshers in general.

---

1. Segment recovery fails on simple polygons with reflex-adjacent corners

import torch
from physicsnemo.mesh import Mesh
from physicsnemo.mesh.tessellation import fill_interior

pts = torch.tensor([[-0.6, -0.3], [-0.4, -0.7], [-0.1, -0.3],
                    [0.1, -1.0], [0.3, -0.7], [0.5, -0.2]], dtype=torch.float64)
i = torch.arange(6)
fill_interior(Mesh(points=pts, cells=torch.stack([i, (i + 1) % 6], 1)))
# RuntimeError: segment recovery could not find a starting wedge; input
# geometry is too degenerate for float64 predicates

Silent form — another exactly-simple hexagon; the picture shows the
triangulation quietly missing 25% of the domain. With min_angle_degrees=30,
this seems to loop indefinitely. With min_angle_degrees=0.0 we get a
poorly formed mesh.

import torch
import matplotlib.pyplot as plt
from physicsnemo.mesh.tessellation.delaunay import delaunay_mesh_2d

pts = torch.tensor([[0.97, 0.00], [0.32, 0.09], [0.43, 0.26],
                    [0.23, 0.44], [0.43, -0.83], [0.84, -0.24]], dtype=torch.float64)
# NOTE: the default min_angle_degrees=30 hangs on this input; 0 isolates the constrained triangulation.
p, t, _, _ = delaunay_mesh_2d([pts], max_area=None, min_angle_degrees=0.0)
x, y = pts[:, 0], pts[:, 1]
exact = float(0.5 * (x * y.roll(-1) - x.roll(-1) * y).sum().abs())
tp = p[t]
meshed = float((0.5 * ((tp[:, 1, 0] - tp[:, 0, 0]) * (tp[:, 2, 1] - tp[:, 0, 1])
       - (tp[:, 2, 0] - tp[:, 0, 0]) * (tp[:, 1, 1] - tp[:, 0, 1]))).abs().sum())
fig, ax = plt.subplots()
ax.fill(x, y, color="tab:red", alpha=0.25, lw=0)  # red = domain the mesh fails to cover
ax.tripcolor(p[:, 0], p[:, 1], t, facecolors=0.0 * t[:, 0], cmap="Greys", vmin=0, vmax=1)
ax.triplot(p[:, 0], p[:, 1], t, lw=0.4, color="tab:blue")
ax.set_aspect("equal")
ax.set_title(f"exactly simple hexagon: area in {exact:.3f}, meshed {meshed:.3f}, no error")
plt.show()

Root cause traced:
_first_crossing picks the starting wedge with two cross-sign tests that
never confirm the ray points into the wedge, so near-straight reflex
corners select the wedge behind the vertex; the pipe walk then exits the hull
and _collect_pipe computes 3*n with n = -1 — Python negative indexing
silently reads garbage adjacency. Suggested guards: pick the wedge by proper
edge-vs-segment intersection, assert n >= 0 in _collect_pipe, cap
_triangles_around. The current test shapes are all densely sampled smooth
boundaries where recovery is a no-op; both hexagons above were found by
shrinking (delta debugging), so a shrinking property-based test over random
star polygons would lock this down.

2. Component grouping misclassifies valid rectilinear nesting (regression in 4c3be42)

import torch
from physicsnemo.mesh import Mesh
from physicsnemo.mesh.tessellation import fill_interior

# 6x6 plate with a centered 5x5 cutout: a valid square annulus.
points = torch.tensor([[0.0, 0.0], [6.0, 0.0], [6.0, 6.0], [0.0, 6.0],    # outer loop
                       [0.5, 0.5], [5.5, 0.5], [5.5, 5.5], [0.5, 5.5]],   # hole loop
                      dtype=torch.float64)
cells = torch.tensor([[0, 1], [1, 2], [2, 3], [3, 0],
                      [4, 5], [5, 6], [6, 7], [7, 4]])
try:
    fill_interior(Mesh(points=points, cells=cells))
except ValueError as e:
    print("annulus crash:", e)   # torch.cat(): expected a non-empty list of Tensors

# Same plate plus a 1x1 island inside the hole: valid domain, area 36 - 25 + 1 = 12.
points = torch.cat([points, torch.tensor([[2.5, 2.5], [3.5, 2.5], [3.5, 3.5], [2.5, 3.5]],
                                          dtype=torch.float64)])
cells = torch.cat([cells, torch.tensor([[8, 9], [9, 10], [10, 11], [11, 8]])])
m = fill_interior(Mesh(points=points, cells=cells))
print(f"island-in-hole: {m.n_points} points, {m.n_cells} cells")
# island-in-hole: 4 points, 2 cells  <- the island alone; 92% of the domain silently dropped

Cause: _group_components probes containment with
polygon_interior_point(loop_j) (largest-ear centroid), and a square outer
loop's ear centroid legitimately falls inside its own hole, so every loop
lands at odd depth. Not a thin-ring corner case — any centered square hole of
side ≥ 2 in the 6x6 plate fails; ≤ 1.8 works. Containment between disjoint
loops should be tested with a point of loop j's boundary (any vertex);
touching loops violate the documented disjointness precondition and can just
raise.

3. Corners sharper than min_angle_degrees silently break the guarantee

import math
import torch
from physicsnemo.mesh.tessellation.delaunay import delaunay_mesh_2d

wedge = torch.tensor([[0., 0.], [1., 0.],
                      [math.cos(math.radians(20)), math.sin(math.radians(20))]],
                     dtype=torch.float64)
p, t, _, _ = delaunay_mesh_2d([wedge], max_area=None, min_angle_degrees=30.0)
tp = p[t]
areas = 0.5 * ((tp[:,1,0]-tp[:,0,0]) * (tp[:,2,1]-tp[:,0,1])
             - (tp[:,2,0]-tp[:,0,0]) * (tp[:,1,1]-tp[:,0,1]))
d = (p - wedge[0]).norm(dim=-1)
print(len(p), "points;", int((d < 1e-9).sum()), "within 1e-9 of the apex")
# 717 points; 597 within 1e-9 of the apex
print("smallest triangle area:", float(areas.abs().min()))
# smallest triangle area: 2.9e-109
t32 = p.float()[t]
a32 = 0.5 * ((t32[:,1,0]-t32[:,0,0]) * (t32[:,2,1]-t32[:,0,1])
           - (t32[:,2,0]-t32[:,0,0]) * (t32[:,1,1]-t32[:,0,1]))
print(int((a32 == 0).sum()), "of", len(t), "triangles have exactly zero area in float32")
# 424 of 715 triangles have exactly zero area in float32

The run terminates "successfully" - no RuntimeError, though the docs promise
one for inputs outside the ≥60° assumption — because _is_bad's raw product
of squared edge lengths underflows to 0 > 0 = False before the insertion
budget trips. A 20° corner makes the 30° bound geometrically unsatisfiable at
the apex, so the honest behaviors are reject-loudly or relax-and-document;
instead the mesh ships with debris that detonates
downstream. Suggested regardless of policy: make _is_bad scale-invariant
so underflow can't mask bad triangles, guard midpoint == endpoint in _split_segment,
and re-scope the ≥60° caveat — it currently reads as applying only to "segments of distinct loops", i.e.
same-loop corners unrestricted, which is exactly the failing case.

[peterdsharpe]

Thanks @melo-gonzo — all three cases were valid, and this was exactly the review the PR needed. Fixed in 3d256a7:

1. Segment recovery (valid — and worse than it looked). Your root-cause diagnosis was correct: the starting-wedge test checked only that the wedge's two vertices straddle the ray line, which can't distinguish the wedge containing the forward ray from the backward one — your second hexagon showed the walk then exits the hull through Python negative indexing and silently drops coverage. Not a general reflex-corner limitation: a simple polygon always admits a CDT, so this was squarely our bug. Wedge selection now requires a proper segment-segment intersection with the wedge's opposite edge (taking the crossing nearest the segment start, since a segment can properly cross a non-convex link polygon more than once), and the pipe walk hard-fails on hull exit. Your observation about why the suite missed it (densely sampled boundaries make recovery a no-op) is now institutionalized: both hexagons are named regression tests, plus a 25-case random-star-polygon property test asserting shoelace-exact area on pure CDTs.

2. Nesting misclassification (valid — a regression from an earlier review round). Exactly as you said: the largest-ear interior point of a square outer loop legitimately falls inside its own hole. Containment between disjoint curves is now probed with the loops' boundary vertices — all of them, so touching/crossing loops raise a mixed-membership error instead of classifying arbitrarily. Your square-annulus and island-in-hole cases are regression tests with exact area assertions.

3. Sharp corners (valid). The ~60° corner condition was documented as an assumption but never enforced; your 20° wedge shows the failure mode is silent and severe (sub-float32-area triangles). It's now validated: a corner sharper than 60° with an angle bound requested raises with the corner's location and the min_angle_degrees=0.0 escape hatch (pure CDT still meshes your wedge exactly). A guaranteed-quality fix for sharp corners (Miller–Pav–Walkington-style concentric shells) is possible future work, but loud rejection beats silent garbage today.

Also scrubbed the notebook execution metadata as suggested.

[peterdsharpe]

/ok to test 3d256a7

[peterdsharpe]

/ok to test d77f87b

[peterdsharpe]

/ok to test 2c15c2b