Geometry

Geometry on a DCEL mesh: trijunction angles and lengths, areas, volumes, and their derivatives.

Computation of trijunction mean angles (in radians, under any of the angle readings of foambryo.second_order_angle_reading) and lengths, interface and region areas, region volumes, and their derivatives with respect to the vertex positions. Every quantity is in the mesh's own length unit; interfaces are keyed (label_1, label_2), label_1 < label_2, label 0 the exterior medium.

Sacha Ichbiah, Sept 2021. Matthieu Perez, 2024.

TrijunctionEdge dataclass

TrijunctionEdge(vertices: tuple[int, int], length: float, half_edge_keys: tuple[int, ...], keys: tuple[tuple[int, int, int], ...], angles: tuple[float, ...])

One geometric trijunctional edge, with the three wedge angles measured on it.

The per-edge counterpart of :func:compute_angles_tri's per-triple-line averages, and the object TensionComputationMethod.YoungDupreLocal assembles its rows from. One instance per geometric edge, not one per half-edge: half_edge_keys records the de-duplicated multiplicity so a caller can check it rather than assume it.

Attributes:

Name	Type	Description
vertices	tuple[int, int]	the edge's two vertex keys, ordered (min, max).
length	float	the edge's length. Not accumulated, not over-counted.
half_edge_keys	tuple[int, ...]	keys of every half-edge incident on this geometric edge, sorted. Its length is the visit multiplicity — 3 for a trijunctional edge on a closed mesh.
keys	tuple[tuple[int, int, int], ...]	the three wedge keys (min(e, g), f, max(e, g)), one per adjacent region f.
angles	tuple[float, ...]	the wedge angle in radians, aligned index-wise with keys.

compute_length_derivative_dict

compute_length_derivative_dict(mesh: DcelData) -> dict[tuple[int, int, int], NDArray[np.float64]]

Compute map trijunction (id reg 1, id reg 2, id reg 3) -> derivative of length wrt each point.

Closed-form float64 replacement for the former autodiff path (moved to foambryo._autodiff_derivatives). For an edge with endpoints p, q the total-length gradient is the unit edge vector at p and its negative at q: d||p - q|| / dp = (p - q) / ||p - q||. Summed over every edge of the trijunction.

compute_length_trijunctions

compute_length_trijunctions(mesh: DcelData, prints: bool = False) -> dict[tuple[int, int, int], float]

Compute the total length of each trijunction of the mesh in a map (id reg 1, id reg 2, id reg 3) -> length.

compute_length_edges_trijunctions

compute_length_edges_trijunctions(points: NDArray[float64], edges_trijunctions: NDArray[int64]) -> NDArray[np.float64]

Compute the length of each edge of a trijunction .

extract_edges_trijunctions

extract_edges_trijunctions(mesh: DcelData, prints: bool = False) -> dict[tuple[int, int, int], NDArray[np.int64]]

Extract a dict that maps trijunctions (id reg 1, id reg 2, id reg 3) to a list of edges.

compute_areas_faces

compute_areas_faces(mesh: DcelData) -> None

Compute area of every triangle in the mesh. Modify the info of the mesh.

compute_areas_cells

compute_areas_cells(mesh: DcelData) -> dict[int, float]

Compute area of each cells.

compute_areas_interfaces

compute_areas_interfaces(mesh: DcelData) -> dict[tuple[int, int], float]

Compute area of every interface (label1, label2) in mesh.

compute_area_derivative_dict

compute_area_derivative_dict(mesh: DcelData) -> dict[tuple[int, int], NDArray[np.float64]]

Compute dict that maps interface (label1, label2) to array of change of area per point.

compute_area_derivative_sparse

compute_area_derivative_sparse(mesh: DcelData) -> tuple[list[tuple[int, int]], csr_matrix]

Sparse area-gradient operator D (n_interfaces x 3*n_points), one row per interface.

Row i is the flattened gradient of interface keys[i]'s total area with respect to every vertex coordinate; vertex v's three components sit at columns 3v, 3v+1, 3v+2. The per-vertex contributions are identical to :func:compute_area_derivative_dict, but the operator is assembled by a single vectorised COO scatter-add over the triangle list (O(n_triangles) non-zeros, mean support ~6 % of points per interface) — no Python loop over interfaces or vertices, and never a dense (n_interfaces, n_points, 3) array. This is the object the variational tension/pressure solve consumes.

Returns:

Type	Description
tuple[list[tuple[int, int]], csr_matrix]	tuple[list[tuple[int, int]], csr_matrix]: the sorted interface keys (row order) and the sparse D.

compute_volume_cells

compute_volume_cells(mesh: DcelData) -> dict[int, float]

Compute map cell number -> volume.

compute_volume_derivative_dict

compute_volume_derivative_dict(mesh: DcelData) -> dict[int, NDArray[np.float64]]

Compute map cell number -> derivative of volume wrt to each point.

compute_volume_derivative_sparse

compute_volume_derivative_sparse(mesh: DcelData) -> tuple[NDArray[np.int64], csr_matrix]

Sparse volume-gradient operator V (n_materials x 3*n_points), one row per material.

Row j is the flattened gradient of material materials[j]'s signed volume. Rows are ordered as mesh.materials (so materials[0] is the exterior); the variational solve uses rows [1:], the cells. Each triangle contributes +cross to its lower-label material and -cross to its higher-label one, exactly as :func:compute_volume_derivative_dict, but assembled by one vectorised COO scatter-add instead of six nested n_materials x n_points Python lists.

Returns:

Type	Description
tuple[NDArray[int64], csr_matrix]	tuple[NDArray[np.int64], csr_matrix]: mesh.materials (row order) and the sparse V.

compute_angles_tri

compute_angles_tri(mesh: DcelData, unique: bool = True, *, reading: str = CHORD_READING) -> tuple[dict[tuple[int, int, int], float], dict[tuple[int, int, int], float], dict[tuple[int, int, int], float]]

Compute three maps trijunction (id reg 1, id reg 2, id reg 3) to mean angle, mean angle (deg), length.

The keys are wedge keys: (min(e, g), f, max(e, g)) is the angle of the wedge inside region f at the triple line where regions e, f and g meet, so one triple line contributes three keys whose angles sum to 360 degrees.

That sum holds exactly only under an oriented reading. Under "chord" and "quadric" the angle is arccos of a dot product, whose range is [0, 180] degrees, so a wedge wider than a straight angle comes back as 360 - theta and its edge's three angles sum to 360 - 2 (theta - 180). "oriented_chord" and "oriented_quadric" return the wedge in [0, 360) and restore the identity; they are bitwise identical to their unoriented counterparts on every wedge that is not reflex.

unique. In earlier versions the parameter had no effect except to double every accumulator uniformly, which changes neither the mean angle nor the length ratio that every builder actually uses. It now means what its name says:

• unique=True (the declared default) counts each geometric trijunctional edge exactly once per wedge key. dict_length[key] is then the true arclength of that triple line, and sum(dict_length.values()) / 3 — the total_length_mesh every builder computes — is then the true total junction length of the mesh, which it was not before.
• unique=False reproduces the historical accumulation bit for bit: each geometric edge is visited once per incident half-edge (3) and each visit appends twice, so dict_length[key] = 6 * length and sum(...) / 3 = 6 * L_total. Every call site inside foambryo.tension_inference passes unique=False, so no inferred number moves.

The mean angles are identical either way — the over-count was, and remains, uniform on the benchmarking dataset (BENCHMARKS.md). That uniformity was verified rather than assumed: on all 47 of its meshes every trijunctional geometric edge has exactly 3 incident half-edges, and the three visits return bit-identical angles. The one construction where that edge accounting is not uniform — a triple line running off an open mesh boundary, where an edge that has lost one of its three films is dropped and the line's length weight is under-stated in proportion — is pinned in tests/test_per_edge_young_dupre.py.

reading. "chord" (the default) is the convention described above; adding the selector changed none of its numbers. "oriented_chord" is the same reading with the wedge angle reported in [0, 360) degrees instead of [0, 180], so a wedge wider than a straight angle is no longer reported as its reflection and the three wedge keys of a line really do sum to 360 degrees; it is bitwise identical to "chord" wherever no wedge is reflex. "quadric" replaces each triangle's chord by the tangent of a generic quadric fitted to that interface's own mesh vertices around the junction — see foambryo.second_order_angle_reading. It is not available on every face of every mesh: where the fit is inadmissible the chord is kept, so a "quadric" result can mix the two conventions. This function reports that mixing as an AngleReadingCoverageWarning; use :func:compute_angles_tri_with_report to receive the coverage numbers instead.

Parameters:

Name	Type	Description	Default
mesh	DcelData	Mesh to analyze.	required
unique	bool	Count each geometric edge once (True) or reproduce the historical 6x over-count (False). Defaults to True.	True
reading	str	"chord", "quadric", "oriented_chord" or "oriented_quadric". Keyword-only. Defaults to "chord", which is bit-identical to the historical behaviour.	CHORD_READING

Returns:

Type	Description
tuple[dict[tuple[int, int, int], float], dict[tuple[int, int, int], float], dict[tuple[int, int, int], float]]	tuple[dict, dict, dict]: mean angle (rad), mean angle (deg) and accumulated length, all keyed by wedge key.

compute_angles_tri_with_report

compute_angles_tri_with_report(mesh: DcelData, unique: bool = True, *, reading: str = CHORD_READING, skip_edges: Collection[tuple[int, int]] | None = None) -> tuple[dict[tuple[int, int, int], float], dict[tuple[int, int, int], float], dict[tuple[int, int, int], float], AngleReadingReport]

:func:compute_angles_tri, plus the reading's coverage report.

The reporting entry point. Use it whenever a result will be published, compared or benchmarked: it is the only way to know which trijunction lines a "quadric" reading could not cover and therefore read as chords.

Parameters:

Name	Type	Description	Default
mesh	DcelData	Mesh to analyze.	required
unique	bool	Count each geometric edge once (True) or reproduce the historical 6x over-count (False). Defaults to True.	True
reading	str	one of foambryo.second_order_angle_reading.ANGLE_READINGS. Keyword-only. Defaults to "chord".	CHORD_READING
skip_edges	Collection[tuple[int, int]] | None	geometric edges, as (min, max) vertex-key pairs, whose wedges contribute neither an angle nor a length to any triple line's accumulators. Keyword-only, None by default, and with None (or an empty collection) this function's result is bitwise what it was before the parameter existed — the membership test is the only thing that runs. Supplied by foambryo.wedge_refusal so that an estimator can decline an observation its own geometry makes infeasible; an edge listed here is simply not measured, so a line that loses every one of its edges disappears from all three returned maps rather than appearing with an empty mean.	None

Returns:

Type	Description
tuple[dict[tuple[int, int, int], float], dict[tuple[int, int, int], float], dict[tuple[int, int, int], float], AngleReadingReport]	tuple[dict, dict, dict, AngleReadingReport]: mean angle (rad), mean angle (deg), accumulated length, and the coverage report.

compute_trijunction_edges

compute_trijunction_edges(mesh: DcelData, *, reading: str = CHORD_READING) -> list[TrijunctionEdge]

List every geometric trijunctional edge with its three wedge angles, each edge once.

The report-dropping wrapper around :func:compute_trijunction_edges_with_report; it emits an AngleReadingCoverageWarning in the report's place when a "quadric" reading fell back to chords anywhere.

Parameters:

Name	Type	Description	Default
mesh	DcelData	Mesh to analyze.	required
reading	str	one of foambryo.second_order_angle_reading.ANGLE_READINGS. Keyword-only. Defaults to "chord".	CHORD_READING

Returns:

Type	Description
list[TrijunctionEdge]	list[TrijunctionEdge]: one record per geometric trijunctional edge, in increasing order of lowest incident half-edge key.

compute_trijunction_edges_with_report

compute_trijunction_edges_with_report(mesh: DcelData, *, reading: str = CHORD_READING) -> tuple[list[TrijunctionEdge], AngleReadingReport]

:func:compute_trijunction_edges, plus the reading's coverage report.

Uses the same convention helper as :func:compute_angles_tri (:func:_wedges_at_half_edge), so averaging this function's angles over the edges of a triple line reproduces compute_angles_tri's mean for that line to machine precision — pinned by tests/test_per_edge_young_dupre.py.

De-duplication is by the geometric edge's (min, max) vertex-key pair, keeping the lowest-keyed incident half-edge; the choice does not matter numerically (the three visits of every trijunctional edge in the 47 meshes of the benchmarking dataset were measured to be bit-identical) but it makes the output deterministic. Quadrijunctional edges are skipped, exactly as compute_angles_tri skips them.

reading selects the angle convention exactly as :func:compute_angles_tri does, and it is the route by which TensionComputationMethod.YoungDupreLocal reaches the second-order reading. The edge set, the lengths and the de-duplication do not depend on it — only the angles do.

Parameters:

Name	Type	Description	Default
mesh	DcelData	Mesh to analyze.	required
reading	str	one of foambryo.second_order_angle_reading.ANGLE_READINGS. Keyword-only. Defaults to "chord".	CHORD_READING

Returns:

Type	Description
tuple[list[TrijunctionEdge], AngleReadingReport]	tuple[list[TrijunctionEdge], AngleReadingReport]: one record per geometric trijunctional edge, in increasing order of lowest incident half-edge key, and which trijunction faces the reading covered.