Validation

Validation and benchmarking tools: the unified force-balance operator and mesh checks, on raw mesh arrays.

The unified force-balance operator (the gradient of the capillary energy with respect to the mesh vertices, one column per unknown), its solvers and diagnostics, so that it can be tested and benchmarked within foambryo; the ForceBalance* tension methods are built on it. The test suite compares these functions with an independent reference implementation of the same algorithms (tests/reference/force_balance_reference_implementation.py). These functions check and benchmark an inference; they are not an inference method of their own.

A mesh here is three arrays, (V, F, regions): vertices (nv, 3), triangles (nf, 3) (vertex indices) and the two regions each triangle separates, (nf, 2), the medium being region 0 by default (foambryo.validation.mesh_io); foambryo.validation.force_balance builds and solves the operator on them.

block_anchors

block_anchors(ni: int, nj: int, ncols: int) -> list[NDArray[np.float64]]

Build the two constraint vectors anchoring the mean surface and line tension.

Passing both to solve_constrained fixes the scale of each block.

Parameters:

Name	Type	Description	Default
ni	int	Number of interfaces (surface-tension columns).	required
nj	int	Number of junctions (line-tension columns).	required
ncols	int	Total number of columns of M.	required

Returns:

Type	Description
list[NDArray[float64]]	list[NDArray[np.float64]]: [a_surface, a_line], each the mean of its block's columns. With nj == 0, a_line is all zeros: pass a_surface alone to solve_constrained then.

constrained_covariance

constrained_covariance(m: NDArray[float64], anchors: list[NDArray[float64]] | NDArray[float64], sigma_f: float) -> NDArray[np.float64]

Covariance of the solve_constrained estimator under i.i.d. force-balance noise.

Model: the vertex force residual M @ beta is not exactly zero but epsilon ~ N(0, sigma_f^2 I), so beta_hat is a constrained generalised-least-squares estimator. Its covariance is the usual GLS covariance projected onto the constraint's null space,

Sigma = sigma_f^2 (H^-1 - H^-1 A^T (A H^-1 A^T)^-1 A H^-1),     H = M^T M,

which is singular along every constraint direction — variance exactly zero in the gauge that was fixed by hand, which is the correct statement: the overall tension scale is a convention, not a measured quantity. This is the error bar that belongs with the KKT solve, and it needs no prior; contrast map_estimate, whose posterior covariance is non-singular only because the prior supplies the missing information.

Honest caveat: geometric noise on a good mesh is correlated between neighbouring vertices, since they share triangles, so an i.i.d. epsilon understates the true uncertainty. Use these as relative, comparable error bars across interfaces, not as calibrated absolute ones.

Parameters:

Name	Type	Description	Default
m	NDArray[float64]	The force-balance operator.	required
anchors	list[NDArray[float64]] | NDArray[float64]	The same constraint row vectors passed to solve_constrained.	required
sigma_f	float	Standard deviation of the per-component vertex force residual. residual_noise_scale estimates it from the fit itself.	required

Returns:

Type	Description
NDArray[float64]	NDArray[np.float64]: Covariance of beta, shape (ncols, ncols).

effective_surface_tensions

effective_surface_tensions(m: NDArray[float64], ni: int, nj: int, mean_tension: float = 1.0) -> tuple[NDArray[np.float64], NDArray[np.float64], float]

Ask which surface tensions balance the mesh, fixing nothing but the overall gauge.

This is the surface-tension analogue of the effective-line-tension check (keep the true surface tensions and ask which line tensions balance the mesh), which detects a factor on a subset of interfaces (for example, some line tensions encoded at half their value by the program that generated the mesh) — a defect in the data, not in the inference. Here the harness operator carries no line-tension block (beta = (gamma, p), by design), so the same question is asked of the surface tensions directly: solve

min ||M beta||^2   subject to   mean(gamma) = mean_tension

and compare gamma to nominal. Only ratios are physical, so the constraint is a gauge choice, not a physical input — without it beta = 0 is the trivial minimiser. A comparison of the returned tensions against nominal ones is therefore only meaningful after gauge matching (see foambryo_benchmarks.metrics.effective_tension_ratios in the source repository).

The only constraint imposed is the mean surface tension, because the harness operator has no line block; with a line block, block_anchors gives a constraint for each tension block. nj > 0 is accepted here only so the residual is reported consistently, and the caller is expected to pass nj = 0.

Parameters:

Name	Type	Description	Default
m	NDArray[float64]	The force-balance operator.	required
ni	int	Number of interfaces (surface-tension columns).	required
nj	int	Number of junctions (line-tension columns).	required
mean_tension	float	The gauge: the mean of the returned surface tensions. Defaults to 1.0, matching foambryo.infer_tensions' own convention and the benchmarking-dataset ground truth.	1.0

Returns:

Type

Description

tuple[NDArray[float64], NDArray[float64], float]

tuple[NDArray[np.float64], NDArray[np.float64], float]: (gamma_eff, p_eff, residual_fraction), where residual_fraction is ||M beta_eff|| / ||M_gamma gamma_eff|| — the force left unbalanced once the tensions themselves are free (rho_free). It is the ratio at the minimiser of ||M beta|| under the gauge, not the minimum of the ratio itself, so it is close to, but not strictly, the lowest value residual_at's rho can take at any given tensions. rho far above it means the nominal tensions are wrong; both being large means the mesh is not an equilibrium of any tensions at all.

effective_tensions

effective_tensions(m: NDArray[float64], ni: int, nj: int, gamma: NDArray[float64]) -> NDArray[np.float64]

Given the surface tensions, find the line tensions (and pressures) that best balance the mesh.

Parameters:

Name	Type	Description	Default
m	NDArray[float64]	The force-balance operator.	required
ni	int	Number of interfaces (surface-tension columns).	required
nj	int	Number of junctions (line-tension columns).	required
gamma	NDArray[float64]	Surface tensions, shape (ni,).	required

Returns:

Type	Description
NDArray[float64]	NDArray[np.float64]: Best-fit line tensions, shape (nj,).

force_balance_operator

force_balance_operator(v: NDArray[float64], f: NDArray[int64], regions: NDArray[int64], interfaces: list[tuple[int, int]], junctions: list[tuple[int, int, int]], edge_to_triple: dict[tuple[int, int], tuple[int, int, int]], cells: list[int]) -> NDArray[np.float64]

Assemble M of shape (3*nv, n_interfaces + n_junctions + n_cells).

Column order: surface tensions, line tensions, pressures. Row i (three components) is the energy gradient at vertex i.

Parameters:

Name	Type	Description	Default
v	NDArray[float64]	Mesh vertices, shape (nv, 3).	required
f	NDArray[int64]	Mesh triangles (vertex indices), shape (nf, 3).	required
regions	NDArray[int64]	Per-triangle region pair (first, second), shape (nf, 2); the triangle's right-hand-rule normal points out of the first region into the second (check_label_orientation_convention).	required
interfaces	list[tuple[int, int]]	Sorted interface keys (a, b), a < b (see list_interfaces).	required
junctions	list[tuple[int, int, int]]	Sorted junction keys (see list_junctions). Pass [] for the no-line-tension operator M = [M_gamma | M_p].	required
edge_to_triple	dict[tuple[int, int], tuple[int, int, int]]	Map mesh edge -> its region-triple (see list_junctions). Pass {} together with junctions=[].	required
cells	list[int]	Sorted cell labels, the medium left out (see cell_labels).	required

Returns:

Type	Description
NDArray[float64]	NDArray[np.float64]: The operator M, shape (3*nv, n_interfaces + n_junctions + n_cells), in the mesh's own units.

Raises:

Type	Description
KeyError	if a triangle's region pair is not in interfaces, or an edge's triple not in junctions.

map_estimate

map_estimate(m: NDArray[float64], prior_mean: NDArray[float64], prior_cov: NDArray[float64], sigma_f: float = 1.0, nonneg_upto: int | None = None) -> tuple[NDArray[np.float64], NDArray[np.float64]]

Gaussian MAP estimate of beta, with posterior covariance.

Minimises |M @ beta|^2 / (2 * sigma_f^2) + 1/2 (beta - prior_mean)^T prior_cov^-1 (beta - prior_mean). A prior with nonzero mean on both tension blocks softly anchors both scales, and the posterior covariance gives error bars.

Parameters:

Name	Type	Description	Default
m	NDArray[float64]	The force-balance operator.	required
prior_mean	NDArray[float64]	Prior mean of beta.	required
prior_cov	NDArray[float64]	Prior covariance of beta.	required
sigma_f	float	Force-balance noise scale. Defaults to 1.0.	1.0
nonneg_upto	int | None	If given (e.g. n_surface + n_line), those first entries are kept non-negative by 500 projected-gradient steps (pressures left free). Defaults to None.	None

Returns:

Type	Description
tuple[NDArray[float64], NDArray[float64]]	tuple[NDArray[np.float64], NDArray[np.float64]]: (beta, posterior_cov); posterior_cov is the Gaussian posterior's, without the non-negativity constraint.

marginalise_pressures

marginalise_pressures(m: NDArray[float64], ni: int, nj: int = 0) -> tuple[NDArray[np.float64], NDArray[np.float64]]

Eliminate the pressure block in closed form, keeping the tension columns.

For any fixed tensions the pressures that minimise ||M beta||^2 are unique (when M_p has full column rank) and given by p* = -M_p^+ (M_gamma gamma + M_t t). Substituting them back leaves

min ||P_perp [M_gamma | M_t] (gamma, t)||^2,     P_perp = I - M_p M_p^+,

an exactly equivalent problem in the tensions alone — this is block elimination of a quadratic form, not an approximation, so any constraint acting only on the tension columns (the mean-tension gauge) gives the same tensions from the reduced and the joint problem. The test test_marginalised_and_joint_solves_agree in tests/test_force_balance_inference.py checks this numerically.

P_perp is applied through a thin QR factorisation of M_p rather than by forming the 3nv x 3nv projector, which would be a ~9 GB dense array on a 20 000-vertex mesh.

Parameters:

Name	Type	Description	Default
m	NDArray[float64]	The force-balance operator.	required
ni	int	Number of interfaces (surface-tension columns).	required
nj	int	Number of junctions (line-tension columns). Defaults to 0, the harness-default operator M = [M_gamma | M_p].	0

Returns:

Type	Description
tuple[NDArray[float64], NDArray[float64]]	tuple[NDArray[np.float64], NDArray[np.float64]]: (reduced, pressure_map). reduced is P_perp [M_gamma | M_t], shape (3nv, ni + nj); pressure_map recovers the eliminated pressures as p = pressure_map @ (gamma, t), shape (n_cells, ni + nj).

residual_at

residual_at(m: NDArray[float64], ni: int, nj: int, gamma: NDArray[float64], t: NDArray[float64]) -> float

Fraction of the surface-tension force left unbalanced once pressures are fitted.

This is the equilibrium residual rho: the norm of the force left over at the given tensions once the best pressures are fitted, relative to the norm of the surface-tension force. With b = M_gamma @ gamma + M_t @ t, the best pressures are p* = -M_p^+ @ b and rho = ||b + M_p @ p*|| / ||M_gamma @ gamma|| = ||P_perp @ b|| / ||M_gamma @ gamma||, with P_perp = I - M_p M_p^+ the projector orthogonal to the pressure column space, computed here via an explicit least-squares fit of the best pressures rather than forming P_perp directly (equivalent, cheaper). rho ~ 0 means the mesh is at equilibrium at the given (gamma, t).

Parameters:

Name	Type	Description	Default
m	NDArray[float64]	The force-balance operator.	required
ni	int	Number of interfaces (surface-tension columns).	required
nj	int	Number of junctions (line-tension columns).	required
gamma	NDArray[float64]	Surface tensions, shape (ni,).	required
t	NDArray[float64]	Line tensions, shape (nj,). Pass np.zeros(0) when nj == 0 (the no-line-tension operator).	required

Returns:

Name	Type	Description
float	float	rho.

residual_force_field

residual_force_field(m: NDArray[float64], ni: int, nj: int, gamma: NDArray[float64], t: NDArray[float64]) -> NDArray[np.float64]

The per-vertex residual force whose norm residual_at reduces to rho.

Same quantity as inside residual_at — P_perp @ b, computed via an explicit least-squares fit of the best pressures — but returned as a force field on the vertices instead of a single norm. residual_at answers "is this mesh at equilibrium"; this answers "and if not, where", which is what distinguishes a mesh whose triple lines are merely discretised from one whose interfaces have not converged.

Parameters:

Name	Type	Description	Default
m	NDArray[float64]	The force-balance operator.	required
ni	int	Number of interfaces (surface-tension columns).	required
nj	int	Number of junctions (line-tension columns).	required
gamma	NDArray[float64]	Surface tensions, shape (ni,).	required
t	NDArray[float64]	Line tensions, shape (nj,). Pass np.zeros(0) when nj == 0.	required

Returns:

Type	Description
NDArray[float64]	NDArray[np.float64]: Residual force per vertex, shape (nv, 3), and the load M_gamma @ gamma it should be compared against is (m[:, :ni] @ gamma).reshape(-1, 3).

residual_noise_scale

residual_noise_scale(m: NDArray[float64], beta: NDArray[float64], n_constraints: int = 1) -> float

Estimate sigma_f from the fit's own residual: ||M beta||^2 / dof.

The number of degrees of freedom is n_rows - (n_unknowns - n_constraints); a constraint gives one unknown back. This is the textbook residual-variance estimator and it removes the only free hyperparameter from the KKT error bars.

Parameters:

Name	Type	Description	Default
m	NDArray[float64]	The force-balance operator.	required
beta	NDArray[float64]	The solution whose residual is used.	required
n_constraints	int	Number of hard constraints imposed. Defaults to 1 (the mean-tension gauge).	1

Returns:

Name	Type	Description
float	float	The estimated per-component force-residual standard deviation.

singular_values

singular_values(m: NDArray[float64]) -> NDArray[np.float64]

Sorted singular values of the column-normalised operator.

Each column is scaled to unit norm first. The smallest value belongs to the overall-scale (gauge) direction.

Parameters:

Name	Type	Description	Default
m	NDArray[float64]	The force-balance operator.	required

Returns:

Type	Description
NDArray[float64]	NDArray[np.float64]: Ascending singular values of M / ||column||.

soft_mode

soft_mode(m: NDArray[float64], ni: int, nj: int) -> tuple[float, NDArray[np.float64]]

Return (sigma_soft, block_fractions) for the worst-determined non-gauge direction.

The direction is the right singular vector of the column-normalised operator with the second-smallest singular value (the smallest being the gauge). sigma_soft is that singular value; block_fractions is how much of the unit vector sits in the (surface, line, pressure) blocks, as the squared norms of its three parts.

Parameters:

Name	Type	Description	Default
m	NDArray[float64]	The force-balance operator.	required
ni	int	Number of interfaces (surface-tension columns).	required
nj	int	Number of junctions (line-tension columns).	required

Returns:

Type	Description
tuple[float, NDArray[float64]]	tuple[float, NDArray[np.float64]]: (sigma_soft, [surface_fraction, line_fraction, pressure_fraction]).

solve_constrained

solve_constrained(m: NDArray[float64], anchors: list[NDArray[float64]] | NDArray[float64], targets: NDArray[float64]) -> NDArray[np.float64]

Minimise |M @ beta|^2 subject to anchors @ beta = targets, via the KKT system.

anchors is a list of constraint vectors (e.g. the mean of a block, see block_anchors). The system solved is [[2 M^T M, A^T], [A, 0]] @ (beta, lambda) = (0, targets), with A the stacked anchors. With both anchors of block_anchors (mean surface and mean line tension), the scale of each tension block is fixed.

Parameters:

Name	Type	Description	Default
m	NDArray[float64]	The force-balance operator.	required
anchors	list[NDArray[float64]] | NDArray[float64]	Constraint row vectors.	required
targets	NDArray[float64]	Right-hand-side values for each constraint.	required

Returns:

Type	Description
NDArray[float64]	NDArray[np.float64]: beta, the constrained least-squares solution.

Raises:

Type	Description
LinAlgError	if the KKT matrix is singular (for example an all-zero anchor, which block_anchors returns as its line anchor when nj == 0).

split

split(m: NDArray[float64], ni: int, nj: int) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]

Split M into its surface, line and pressure column blocks.

Parameters:

Name	Type	Description	Default
m	NDArray[float64]	The force-balance operator.	required
ni	int	Number of interfaces (surface-tension columns).	required
nj	int	Number of junctions (line-tension columns).	required

Returns:

Type	Description
tuple[NDArray[float64], NDArray[float64], NDArray[float64]]	tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]: (M_gamma, M_t, M_p).

cell_labels

cell_labels(regions: NDArray[int64], exterior: int = 0) -> list[int]

Sorted list of cell regions (everything except the external medium).

Parameters:

Name	Type	Description	Default
regions	NDArray[int64]	Per-triangle region pair, shape (nf, 2).	required
exterior	int	Label of the external medium. Defaults to 0.	0

Returns:

Type	Description
list[int]	list[int]: Sorted cell labels, excluding exterior.

check_label_orientation_convention

check_label_orientation_convention(v: NDArray[float64], f: NDArray[int64], regions: NDArray[int64], exterior: int = 0) -> dict

Check that a mesh's triangle winding agrees with its label column order.

foambryo requires each triangle's right-hand-rule normal to point away from its first-listed region and into its second (see foambryo.dcel._normalize_label_order). That tie between winding and label order cannot be checked from the labels — sorted label columns are neither necessary nor sufficient — so this checks its geometric consequence instead: every bounded cell must enclose a positive volume.

Use this once on any mesh from a source other than dw3d (a simulation, a hand-written .rec, or a mesh whose label columns were sorted by a tool that did not reorient the triangles). Without it the failure is silent: compute_volume_cells simply returns negative volumes and everything downstream consumes them.

Assumes a closed mesh. An open or cropped mesh will report non-positive volumes for the cells it truncates, which is a limitation of the check and not a defect in the mesh; n_nonpositive_cells and the returned volumes let a caller judge which case it is looking at.

Parameters:

Name	Type	Description	Default
v	NDArray[float64]	Mesh vertices, shape (nv, 3).	required
f	NDArray[int64]	Mesh triangles (vertex indices), shape (nf, 3).	required
regions	NDArray[int64]	Per-triangle region pair, shape (nf, 2).	required
exterior	int	Label of the external medium, excluded from the positivity requirement since its volume is legitimately negative. Defaults to 0.	0

Returns:

Name	Type	Description
dict	dict	{"ok", "volumes", "nonpositive_cells", "n_nonpositive_cells", "n_rows_label_unsorted", "n_rows_same_label_both_sides"}. ok is True when every bounded cell has a strictly positive volume.

is_intact

is_intact(v: NDArray[float64], f: NDArray[int64]) -> bool

Check that every face points to an existing vertex (no dangling reference).

Parameters:

Name	Type	Description	Default
v	NDArray[float64]	Mesh vertices, shape (nv, 3).	required
f	NDArray[int64]	Mesh triangles (vertex indices), shape (nf, 3).	required

Returns:

Name	Type	Description
bool	bool	True if no vertex index in f is out of range for v.

list_interfaces

list_interfaces(regions: NDArray[int64]) -> list[tuple[int, int]]

List the region pairs present in the mesh; one surface-tension unknown each.

Parameters:

Name	Type	Description	Default
regions	NDArray[int64]	Per-triangle region pair, shape (nf, 2).	required

Returns:

Type	Description
list[tuple[int, int]]	list[tuple[int, int]]: Sorted, sorted-pair interface keys (a, b) with a < b.

list_junctions

list_junctions(f: NDArray[int64], regions: NDArray[int64]) -> tuple[list[tuple[int, int, int]], dict[tuple[int, int], tuple[int, int, int]]]

Find triple lines: edges carrying three interfaces on exactly three incident triangles.

An edge is admitted as a genuine triple line only under the rule mesh_integrity already documents and enforces as triple_line_ok: three distinct region labels (regions_per_edge == 3), spanning three distinct interfaces (interfaces_per_edge == 3), on exactly three incident triangles (triangles_per_edge == 3) -- not merely three region labels. An edge with three region labels but a fourth (or more) incident triangle is a malformed triple line (mesh_integrity's n_malformed_triple_line_edges), not a genuine one, and is excluded: admitted, it would give a region-triple whose derived interface pair need not be any triangle's actual region pair in the mesh, which junction_graph.partition_by_junction_graph's union-find has no key for (a KeyError).

Parameters:

Name	Type	Description	Default
f	NDArray[int64]	Mesh triangles (vertex indices), shape (nf, 3).	required
regions	NDArray[int64]	Per-triangle region pair, shape (nf, 2).	required

Returns:

Type	Description
tuple[list[tuple[int, int, int]], dict[tuple[int, int], tuple[int, int, int]]]	tuple[list[tuple[int, int, int]], dict[tuple[int, int], tuple[int, int, int]]]: - sorted region-triples, one per distinct genuine triple line; - map from mesh edge (u, v) (u < v) to its region-triple, genuine triple lines only.

mean_edge_length

mean_edge_length(v: NDArray[float64], f: NDArray[int64]) -> float

Average triangle-edge length; handy to express noise or tolerances.

Parameters:

Name	Type	Description	Default
v	NDArray[float64]	Mesh vertices, shape (nv, 3).	required
f	NDArray[int64]	Mesh triangles (vertex indices), shape (nf, 3).	required

Returns:

Name	Type	Description
float	float	Mean length of all (non-deduplicated) triangle edges.

mesh_integrity

mesh_integrity(v: NDArray[float64], f: NDArray[int64], regions: NDArray[int64]) -> dict

Structural integrity report for a multimaterial mesh.

Implements the dataset-integrity check — no dangling vertex references, and every triple line shared by exactly 3 interfaces — and reports the surrounding edge statistics that say how a mesh fails when it does.

intact is deliberately narrow: exactly those two conditions. Everything else (open boundary edges, quadrijunction edges, non-manifold interface edges, degenerate or duplicated triangles) is reported separately, because those are mesh-quality facts rather than integrity failures. In particular quadrijunction edges are expected meshing artifacts, not corruption (in a dry foam exactly three films meet along an edge, so a valence-4 edge is an artifact of discretisation).

Edge terminology, all keyed by the number of distinct region labels appearing on the triangles incident to an edge. This region-label count alone is only a classification, not a validity check: list_junctions uses it to name candidate triple-line edges but, like this function's own n_bad_triple_lines below, only admits one once it also carries three interfaces on exactly three incident triangles (triple_line_ok). A non-manifold edge with three region labels but a fourth incident triangle is malformed, and is left out of the tension system.

• 1 region: an edge interior to a single interface's triangulation;
• 2 regions: an interface edge;
• 3 regions: a triple-line edge — the physical Plateau junction;

• = 4 regions: a quadrijunction edge.

A triple-line edge spanning regions {a, b, c} must carry all three interfaces (a,b), (a,c), (b,c), one triangle each. Fewer means the junction is malformed: the operator would then assemble a force balance for a junction that is not one.

Parameters:

Name	Type	Description	Default
v	NDArray[float64]	Mesh vertices, shape (nv, 3).	required
f	NDArray[int64]	Mesh triangles (vertex indices), shape (nf, 3).	required
regions	NDArray[int64]	Per-triangle region pair, shape (nf, 2).	required

Returns:

Name	Type	Description
dict	dict	Integrity report. intact is True when there is no dangling vertex reference and every triple-line edge carries exactly its three interfaces.

read_arec

read_arec(path: str | Path) -> tuple[NDArray[np.float64], NDArray[np.int64], NDArray[np.int64]]

Read a text multimaterial mesh (.arec).

Parameters:

Name	Type	Description	Default
path	str | Path	Path to the .arec file.	required

Returns:

Type	Description
tuple[NDArray[float64], NDArray[int64], NDArray[int64]]	tuple[NDArray[np.float64], NDArray[np.int64], NDArray[np.int64]]: (V, F, regions).

read_mesh

read_mesh(path: str | Path) -> tuple[NDArray[np.float64], NDArray[np.int64], NDArray[np.int64]]

Read a mesh, choosing the reader from the extension.

Falls back to a binary/text sniff of the first bytes when the extension is neither .rec nor .arec. A .rec is read by read_rec, text or binary.

Parameters:

Name	Type	Description	Default
path	str | Path	Path to the mesh file.	required

Returns:

Type	Description
tuple[NDArray[float64], NDArray[int64], NDArray[int64]]	tuple[NDArray[np.float64], NDArray[np.int64], NDArray[np.int64]]: (V, F, regions).

read_rec

read_rec(path: str | Path) -> tuple[NDArray[np.float64], NDArray[np.int64], NDArray[np.int64]]

Read a multimaterial mesh (.rec), text or binary (told apart by _rec_is_text).

A text .rec (delaunay-watershed-3d's default, and what foambryo mesh writes) is read as read_arec reads an .arec; a binary one in the layout of the module docstring. Before 0.7.0 only the binary layout was read.

Parameters:

Name	Type	Description	Default
path	str | Path	Path to the .rec file.	required

Returns:

Type	Description
tuple[NDArray[float64], NDArray[int64], NDArray[int64]]	tuple[NDArray[np.float64], NDArray[np.int64], NDArray[np.int64]]: (V, F, regions).

require_label_orientation_convention

require_label_orientation_convention(v: NDArray[float64], f: NDArray[int64], regions: NDArray[int64], exterior: int = 0) -> None

Raise if a mesh's winding does not agree with its label column order.

The raising form of check_label_orientation_convention; see that function for what is checked, why label order alone cannot settle it, and the closed-mesh assumption.

Parameters:

Name	Type	Description	Default
v	NDArray[float64]	Mesh vertices, shape (nv, 3).	required
f	NDArray[int64]	Mesh triangles (vertex indices), shape (nf, 3).	required
regions	NDArray[int64]	Per-triangle region pair, shape (nf, 2).	required
exterior	int	Label of the external medium. Defaults to 0.	0

Raises:

Type	Description
ValueError	If any bounded cell encloses a non-positive volume.

signed_cell_volumes

signed_cell_volumes(v: NDArray[float64], f: NDArray[int64], regions: NDArray[int64]) -> dict[int, float]

Divergence-theorem volume per region, with foambryo's sign convention.

Accumulates +det([v0, v1, v2]) / 6 into each triangle's first-listed region and -det(...) / 6 into its second — the same rule as foambryo.geometry.compute_volume_cells, applied to raw arrays so it can be run on a mesh file before building a DcelData.

The result is invariant under DcelData's label-order normalisation (swapping a row's two labels also reverses its winding, which negates det, so the contribution is unchanged), and it was verified equal to DcelData.compute_volumes_cells() on benchmarking-dataset/000_mesh.rec to floating-point precision.

A correctly-wound closed mesh gives a positive volume for every bounded cell and a negative one for the unbounded exterior label, whose faces are all wound inward from its point of view.

Parameters:

Name	Type	Description	Default
v	NDArray[float64]	Mesh vertices, shape (nv, 3).	required
f	NDArray[int64]	Mesh triangles (vertex indices), shape (nf, 3).	required
regions	NDArray[int64]	Per-triangle region pair, shape (nf, 2).	required

Returns:

Type	Description
dict[int, float]	dict[int, float]: Region label -> signed volume.