Tessellation#

Decompose non-simplicial cells into the simplices that the Mesh data structure stores. Currently this provides triangulation of a polygon soup (triangulate()): a vectorized vertex-0 fan for convex polygons and ear clipping for the rare non-convex ones.

A polygon soup is a flat list of polygonal cells, each given only by its vertex ring, with no shared-edge connectivity between cells – the form a surface mesh takes when read straight from a polygonal face list (e.g. a VTP or STL file). Equivalently, it is a surface mesh whose cells may be arbitrary polygons rather than triangles.

Note

This is unrelated to the (now deprecated) Tessellation geometry class in PhysicsNeMo-Sym (physicsnemo.sym.geometry.tessellation), which represents a tessellated STL surface as a solid for sampling surface/interior points and signed-distance values. The utilities here instead split the polygonal cells of an existing mesh into simplices.

Handling non-convex polygons correctly matters for any unsigned-area-weighted quantity (wall-shear / viscous force integration, or total wetted area): the signed vector area of a vertex-0 fan telescopes to the polygon’s total area regardless of convexity, but the sum of the unsigned triangle areas does not.

Every k-gon yields exactly k - 2 triangles, so per-polygon data is broadcast to the output identically in both paths using the returned parent_index.

import torch
from physicsnemo.mesh import Mesh
from physicsnemo.mesh.neighbors import Adjacency
from physicsnemo.mesh.tessellation import triangulate

# A polygon soup as a cell-to-vertex Adjacency (CSR): one quad (vertices 0-3).
points = torch.tensor([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0],
                       [1.0, 1.0, 0.0], [0.0, 1.0, 0.0]])
polygons = Adjacency(offsets=torch.tensor([0, 4]),
                     indices=torch.tensor([0, 1, 2, 3]))

# Low-level: triangle connectivity plus the parent-polygon index.
cells, parent_index = triangulate(points, polygons)
cells          # tensor([[0, 1, 2], [0, 2, 3]])
parent_index   # tensor([0, 0]); broadcast data via cell_data[parent_index]

# High-level: build a Mesh directly, broadcasting per-polygon cell data.
mesh = Mesh.from_polygons(
    points, polygons, cell_data={"pressure": torch.tensor([2.5])}
)

A Mesh can also be constructed in one step with from_polygons().

Exact-Boundary Interior Filling#

Beyond decomposing existing cells, the package also generates quality meshes. fill_interior() takes a closed codimension-one boundary Mesh and fills the enclosed interior with quality simplices. In 2D, the boundary is an edge mesh forming one or more loops, in any order and orientation. The algorithm resolves holes, multiple components, and islands-inside-holes automatically through containment.

The algorithm uses constrained Delaunay triangulation (Bowyer–Watson insertion with constrained-edge recovery, holes removed topologically by even-odd parity flood fill) followed by Ruppert’s Delaunay refinement, so that every 2D output triangle provably satisfies the requested minimum-angle bound and, optionally, a maximum cell size.

Optional optimal-Delaunay-triangulation (ODT) smoothing, controlled by smooth_iterations, moves each interior vertex to the area-weighted average of its incident triangles’ circumcenters (Chen and Xu 2004). This improves the typical angle while preserving both bounds. The exact-boundary contract provides three guarantees:

  • Every input vertex appears bit-identically in the output, in leading rows and input order.

  • The pipeline only ever subdivides boundary facets and never moves them.

  • The entire process is deterministic.

You can attach provenance fields ("boundary_marker" and "source_point") to the output’s point_data by setting provenance=True. By default, the pipeline claims no keys in the user-owned namespace.

The contract is dimension-generic by design. The n = 3 case (a watertight surface Mesh[2, 3] producing tetrahedra) currently raises NotImplementedError: exact 3D boundary recovery is a substantially harder problem that requires its own implementation effort. For implicit domains, or approximate volume meshing of a surface through its SDF, refer to physicsnemo.mesh.generate.mesh_implicit_domain().

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

square = torch.tensor([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
edges = torch.tensor([[0, 1], [1, 2], [2, 3], [3, 0]])
h = 0.1  # target interior edge length
filled = fill_interior(
    Mesh(points=square, cells=edges),
    max_cell_size=math.sqrt(3.0) / 4.0 * h * h,
    min_angle_degrees=30.0,
    smooth_iterations=3,  # optional ODT smoothing
)

polygon_interior_point() is a companion utility returning a point strictly inside a simple polygon.

API Reference#

Tessellation: decompose non-simplicial cells into a simplex mesh.

Exposes triangulate(), which converts a polygon soup (an Adjacency of cell-to-vertex rings) into a triangle connectivity, staying correct for non-convex polygons. It branches on manifold dimension and currently implements the 2D (polygon -> triangle) case.

Also exposes exact-boundary quality mesh generation (which inserts new points, unlike triangulate()’s pure decomposition): fill_interior() takes a closed codimension-one boundary Mesh (edge loops in 2D) and fills the enclosed interior with quality simplices through constrained Delaunay triangulation with Ruppert refinement. The algorithm preserves every input vertex bit-identically, and in 2D every output angle meets the requested bound. polygon_interior_point() returns a point strictly inside a simple polygon.

physicsnemo.mesh.tessellation.fill_interior(
boundary: Mesh,
*,
max_cell_size: float | None = None,
min_angle_degrees: float = 30.0,
smooth_iterations: int = 0,
provenance: bool = False,
) → Mesh[source]#

Fill the interior of a closed boundary mesh with quality simplices.

Dimension-generic contract: given a closed codimension-one boundary Mesh[n-1, n], produce a volume Mesh[n, n] of the enclosed interior such that

  • every input vertex referenced by an edge appears bit-identically in the output (the leading rows of points, in input order);

  • boundary facets are never moved off the input geometry — refinement may subdivide them, but the union of output boundary facets equals the input boundary exactly;

  • interior (Steiner) vertices are inserted to meet the quality bounds.

Currently implemented for n = 2 (edge loops -> triangles), where every output triangle is guaranteed a minimum angle of min_angle_degrees (Ruppert refinement; deterministic, bitwise reproducible). n = 3 (watertight surface -> tetrahedra) raises NotImplementedError pending exact 3D boundary recovery.

Parameters:

  • boundary (Mesh) – Closed codimension-one boundary: n_manifold_dims == n_spatial_dims - 1. In 2D, an edge mesh forming one or more disjoint simple loops (any orientation, any order); nesting is resolved automatically — loops at even containment depth bound components, loops one level deeper bound holes, islands inside holes are supported. Vertices not referenced by any edge are ignored. Loops of distinct components/holes must be disjoint, and segments of distinct loops should meet at angles of at least ~60° for the refinement termination guarantee.

  • max_cell_size (float, optional) – Maximum cell measure (area in 2D). None disables the size bound. For a target edge length \(h\), pass the equilateral measure \(\sqrt{3}/4\,h^2\).

  • min_angle_degrees (float, default 30.0) – Guaranteed minimum triangle angle, in \([0, 33]\) (2D).

  • smooth_iterations (int, default 0) – Quality-gated optimal-Delaunay-triangulation (ODT) smoothing passes after refinement (each interior vertex moves to the area-weighted average of its incident triangles’ circumcenters, Chen & Xu 2004); boundary vertices never move, and the quality bounds are preserved.

  • provenance (bool, default False) –

    When True, attach provenance fields to the output’s point_data (opt-in, since these claim keys in a user-owned namespace):

    • "boundary_marker" (int64): 1 for vertices on the input boundary (input vertices and refinement midpoints inserted on it), 0 for interior Steiner vertices.

    • "source_point" (int64): for vertices inherited from the input, the index of the originating input vertex; -1 for generated vertices. Use it to propagate input point_data onto the output.

Returns:

Volume mesh on the input’s device and dtype, positively oriented, with an empty point_data unless provenance=True.

Return type:

Mesh

Raises:

  • ValueError – If the boundary is not a closed 1-manifold, loops are degenerate or crossing, or quality parameters are out of range.

  • NotImplementedError – For n_spatial_dims != 2 (3D tetrahedralization is planned; the contract above is dimension-generic by design).

Examples

Fill an annulus given as one edge mesh containing both circles:

>>> import math, torch
>>> from physicsnemo.mesh import Mesh
>>> from physicsnemo.mesh.tessellation import fill_interior
>>> def circle(r, n, start):
...     t = torch.arange(n, dtype=torch.float64) / n * 2 * math.pi
...     pts = torch.stack([r * torch.cos(t), r * torch.sin(t)], dim=1)
...     e = torch.stack([torch.arange(n), (torch.arange(n) + 1) % n], dim=1)
...     return pts, e + start
>>> p1, e1 = circle(1.0, 32, 0)
>>> p2, e2 = circle(0.4, 16, 32)
>>> ring = Mesh(points=torch.cat([p1, p2]), cells=torch.cat([e1, e2]))
>>> filled = fill_interior(ring, max_cell_size=0.02)
>>> (filled.n_manifold_dims, filled.n_spatial_dims)
(2, 2)
>>> bool(torch.equal(filled.points[:48], ring.points))  # exact boundary
True
physicsnemo.mesh.tessellation.polygon_interior_point(
loop: Float[Tensor, 'n 2'] | ndarray,
) → Float[Tensor, '2'][source]#

Return a point strictly inside a simple closed polygon.

Returns the centroid of the polygon’s largest ear: a convex vertex whose neighbor-to-neighbor triangle strictly contains no other polygon vertex, which the two-ears theorem guarantees to exist for every simple polygon. Ear triangles lie inside the polygon, so their centroids do too, and taking the largest keeps the point robustly away from the boundary even for polygons with near-degenerate (collinear) vertices. As a final guard against exactly-degenerate inputs, the centroid is verified by an even-odd ray-crossing test before being returned, falling through to the next-largest ear on failure.

The scan is vectorized NumPy: O(n) memory and O(n^2) worst-case time, near-linear in practice because candidate ears are tested largest-first. Deterministic: identical inputs return bitwise-identical outputs.

Parameters:

loop (torch.Tensor or numpy.ndarray) – Polygon vertices, shape \((N, 2)\) with \(N \geq 3\), closed implicitly (do not repeat the first vertex). Must be simple; either winding is accepted.

Returns:

A point strictly inside the polygon, shape \((2,)\), float64.

Return type:

torch.Tensor

Raises:

ValueError – If the loop has fewer than 3 vertices, a non-(N, 2) shape, non-finite coordinates, duplicate consecutive vertices, zero area (fully degenerate), or no verifiable ear (the boundary self-intersects).

physicsnemo.mesh.tessellation.triangulate(
points: Float[Tensor, 'n_points n_spatial'],
polygons: Adjacency,
*,
manifold_dim: int = 2,
assume_convex: bool = False,
) → tuple[Int[Tensor, 'n_simplices d_plus_one'], Int[Tensor, 'n_simplices']][source]#

Decompose cells into simplices, branching on manifold dimension.

Parameters:

  • points (torch.Tensor) – Vertex coordinates, shape \((N_\text{points}, D)\) with \(D \in \{2, 3\}\).

  • polygons (Adjacency) – Cell-to-vertex incidence in CSR form: cell c is the vertex ring polygons.indices[polygons.offsets[c] : polygons.offsets[c + 1]]. Build one from a flat VTK-style soup with Adjacency(offsets=..., indices=connectivity).

  • manifold_dim (int, default 2) – Dimension of the cells to decompose. Only 2 (polygon -> triangle) is implemented.

  • assume_convex (bool, default False) – If True, skip the convexity test and ear-clip fallback and fan-triangulate every cell. Correct only when all cells are convex; this is the fully torch.compile-traceable fast path.

Returns:

  • cells (torch.Tensor) – Simplex connectivity, shape \((N_\text{simplices}, \text{manifold\_dim} + 1)\), dtype int64.

  • parent_index (torch.Tensor) – Source cell of each simplex, shape \((N_\text{simplices},)\). Broadcast per-cell data to the simplices with data[parent_index].

Raises:

  • NotImplementedError – If manifold_dim != 2.

  • ValueError – If any polygon has fewer than three vertices, or if polygons.indices contains a negative index or one >= n_points (i.e. outside the valid \([0, N_\text{points})\) range; both checked off the torch.compile path).

Notes

Each polygon ring must be a simple polygon (no self-intersections) with no repeated consecutive vertices or zero-length edges. Degenerate or self-intersecting rings are not detected and produce undefined results. The non-convex (ear-clip) path additionally assumes each ring is approximately planar: a badly non-planar non-convex ring can project to a self-intersecting 2-D polygon and triangulate incorrectly (convex rings are unaffected by planarity).

Examples

>>> import torch
>>> from physicsnemo.mesh.neighbors import Adjacency
>>> from physicsnemo.mesh.tessellation import triangulate
>>> points = torch.tensor([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0],
...                        [1.0, 1.0, 0.0], [0.0, 1.0, 0.0]])
>>> polygons = Adjacency(offsets=torch.tensor([0, 4]),  # one quad
...                      indices=torch.tensor([0, 1, 2, 3]))
>>> cells, parent_index = triangulate(points, polygons)
>>> cells.tolist()
[[0, 1, 2], [0, 2, 3]]
>>> parent_index.tolist()
[0, 0]