#![forbid(unsafe_code)]
use crate::core::algorithms::incremental_insertion::{
DelaunayRepairErrorKind, InsertionError, InsertionErrorSourceKind,
};
use crate::core::cell::Cell;
use crate::core::collections::{FastHashMap, PeriodicOffsetBuffer, Uuid, VertexKeySet};
use crate::core::operations::InsertionOutcome;
use crate::core::tds::{
CellKey, InvariantError, InvariantErrorSummaryDetail, Tds, TdsError, TdsErrorKind,
TriangulationConstructionState, TriangulationValidationErrorKind, VertexKey,
};
use crate::core::traits::data_type::DataType;
use crate::core::triangulation::{TopologyGuarantee, TriangulationConstructionError};
use crate::core::util::periodic_facet_key_from_lifted_vertices;
use crate::core::vertex::Vertex;
use crate::geometry::kernel::{AdaptiveKernel, Kernel};
use crate::geometry::point::Point;
use crate::geometry::traits::coordinate::{Coordinate, CoordinateScalar};
use crate::topology::spaces::toroidal::ToroidalSpace;
use crate::topology::traits::global_topology_model::{
GlobalTopologyModel, GlobalTopologyModelError,
};
use crate::topology::traits::topological_space::{
GlobalTopology, TopologyKind, ToroidalConstructionMode,
};
use crate::triangulation::delaunay::{
ConstructionOptions, DelaunayRepairPolicy, DelaunayTriangulation,
DelaunayTriangulationConstructionError, DelaunayTriangulationValidationError,
InitialSimplexStrategy, RetryPolicy,
};
use num_traits::ToPrimitive;
use rand::SeedableRng;
use rand::rngs::StdRng;
use rand::seq::SliceRandom;
use std::num::NonZeroUsize;
use thiserror::Error;
const TWO_POW_52_I64: i64 = 4_503_599_627_370_496; const TWO_POW_52_F64: f64 = 4_503_599_627_370_496.0; const MAX_OFFSET_UNITS: i64 = 1_048_576;
const IMAGE_JITTER_UNITS: i64 = 64;
const FNV_OFFSET_BASIS: u64 = 0xcbf2_9ce4_8422_2325;
const FNV_PRIME: u64 = 0x0100_0000_01b3;
type LiftedVertex<const D: usize> = (VertexKey, [i8; D]);
type SymbolicSignature<const D: usize> = Vec<LiftedVertex<D>>;
type PeriodicFacetKey = u64;
type PeriodicCandidate<const D: usize> = (
SymbolicSignature<D>,
SymbolicSignature<D>,
Vec<PeriodicFacetKey>,
bool,
);
fn sort_candidates_by_rarity_and_domain(
order: &mut [usize],
candidate_edges: &[[usize; 3]],
candidate_in_domain: &[bool],
edge_count: usize,
) {
let mut edge_frequency = vec![0usize; edge_count];
for edges in candidate_edges {
for &edge in edges {
edge_frequency[edge] = edge_frequency[edge].saturating_add(1);
}
}
order.sort_by(|a, b| {
let a_edges = candidate_edges[*a];
let b_edges = candidate_edges[*b];
let a_score =
edge_frequency[a_edges[0]] + edge_frequency[a_edges[1]] + edge_frequency[a_edges[2]];
let b_score =
edge_frequency[b_edges[0]] + edge_frequency[b_edges[1]] + edge_frequency[b_edges[2]];
candidate_in_domain[*b]
.cmp(&candidate_in_domain[*a])
.then_with(|| a_score.cmp(&b_score))
.then_with(|| a.cmp(b))
});
}
struct ClosedSelectionDfs<'a> {
target_faces: usize,
order: &'a [usize],
candidate_edges: &'a [[usize; 3]],
edge_counts: Vec<u8>,
selected: Vec<bool>,
nodes: usize,
node_limit: usize,
}
impl ClosedSelectionDfs<'_> {
fn search(&mut self, pos: usize, chosen: usize) -> bool {
if chosen == self.target_faces {
return true;
}
if pos == self.order.len() {
return false;
}
if chosen + (self.order.len() - pos) < self.target_faces {
return false;
}
if self.nodes >= self.node_limit {
return false;
}
self.nodes = self.nodes.saturating_add(1);
let remaining_capacity: usize = self
.edge_counts
.iter()
.map(|&count| usize::from(2_u8.saturating_sub(count)))
.sum();
if chosen + (remaining_capacity / 3) < self.target_faces {
return false;
}
let idx = self.order[pos];
let edges = self.candidate_edges[idx];
if self.edge_counts[edges[0]] < 2
&& self.edge_counts[edges[1]] < 2
&& self.edge_counts[edges[2]] < 2
{
self.selected[idx] = true;
self.edge_counts[edges[0]] += 1;
self.edge_counts[edges[1]] += 1;
self.edge_counts[edges[2]] += 1;
if self.search(pos + 1, chosen + 1) {
return true;
}
self.edge_counts[edges[0]] -= 1;
self.edge_counts[edges[1]] -= 1;
self.edge_counts[edges[2]] -= 1;
self.selected[idx] = false;
}
self.search(pos + 1, chosen)
}
}
fn search_closed_2d_selection(
candidate_edges: &[[usize; 3]],
candidate_in_domain: &[bool],
target_faces: usize,
edge_count: usize,
node_limit: usize,
) -> Option<Vec<bool>> {
let m = candidate_edges.len();
if m < target_faces {
return None;
}
let mut order: Vec<usize> = (0..m).collect();
sort_candidates_by_rarity_and_domain(
&mut order,
candidate_edges,
candidate_in_domain,
edge_count,
);
let mut dfs = ClosedSelectionDfs {
target_faces,
order: &order,
candidate_edges,
edge_counts: vec![0u8; edge_count],
selected: vec![false; m],
nodes: 0,
node_limit,
};
dfs.search(0, 0).then_some(dfs.selected)
}
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
#[non_exhaustive]
pub enum ExplicitTdsErrorKind {
InvalidVertex,
InvalidCell,
InvalidNeighbors,
OrientationViolation,
DuplicateCells,
FacetSharingViolation,
FailedToCreateCell,
NotNeighbors,
MappingInconsistency,
VertexKeyRetrievalFailed,
CellNotFound,
VertexNotFound,
DimensionMismatch,
IndexOutOfBounds,
InconsistentDataStructure,
Geometric,
FacetError,
DuplicateCoordinatesInCell,
}
#[must_use]
#[derive(Clone, Debug, Error, PartialEq, Eq)]
#[error("{message}")]
pub struct ExplicitTdsError {
pub kind: ExplicitTdsErrorKind,
pub message: String,
}
impl From<TdsError> for ExplicitTdsError {
fn from(source: TdsError) -> Self {
let kind = match &source {
TdsError::InvalidVertex { .. } => ExplicitTdsErrorKind::InvalidVertex,
TdsError::InvalidCell { .. } => ExplicitTdsErrorKind::InvalidCell,
TdsError::InvalidNeighbors { .. } => ExplicitTdsErrorKind::InvalidNeighbors,
TdsError::OrientationViolation { .. } => ExplicitTdsErrorKind::OrientationViolation,
TdsError::DuplicateCells { .. } => ExplicitTdsErrorKind::DuplicateCells,
TdsError::FacetSharingViolation { .. } => ExplicitTdsErrorKind::FacetSharingViolation,
TdsError::FailedToCreateCell { .. } => ExplicitTdsErrorKind::FailedToCreateCell,
TdsError::NotNeighbors { .. } => ExplicitTdsErrorKind::NotNeighbors,
TdsError::MappingInconsistency { .. } => ExplicitTdsErrorKind::MappingInconsistency,
TdsError::VertexKeyRetrievalFailed { .. } => {
ExplicitTdsErrorKind::VertexKeyRetrievalFailed
}
TdsError::CellNotFound { .. } => ExplicitTdsErrorKind::CellNotFound,
TdsError::VertexNotFound { .. } => ExplicitTdsErrorKind::VertexNotFound,
TdsError::DimensionMismatch { .. } => ExplicitTdsErrorKind::DimensionMismatch,
TdsError::IndexOutOfBounds { .. } => ExplicitTdsErrorKind::IndexOutOfBounds,
TdsError::InconsistentDataStructure { .. } => {
ExplicitTdsErrorKind::InconsistentDataStructure
}
TdsError::Geometric(_) => ExplicitTdsErrorKind::Geometric,
TdsError::FacetError(_) => ExplicitTdsErrorKind::FacetError,
TdsError::DuplicateCoordinatesInCell { .. } => {
ExplicitTdsErrorKind::DuplicateCoordinatesInCell
}
};
Self {
kind,
message: source.to_string(),
}
}
}
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
#[non_exhaustive]
pub enum ExplicitInsertionErrorKind {
ConflictRegion,
Location,
CavityFilling,
NeighborWiring,
NonManifoldTopology,
HullExtension,
DelaunayValidationFailed,
DelaunayRepairFailed,
DuplicateCoordinates,
DuplicateUuid,
TopologyValidation,
TopologyValidationFailed,
}
#[must_use]
#[derive(Clone, Debug, Error, PartialEq, Eq)]
#[error("{message}")]
pub struct ExplicitInsertionError {
pub kind: ExplicitInsertionErrorKind,
pub source_kind: Option<InsertionErrorSourceKind>,
pub message: String,
}
impl From<InsertionError> for ExplicitInsertionError {
fn from(source: InsertionError) -> Self {
let kind = match &source {
InsertionError::ConflictRegion(_) => ExplicitInsertionErrorKind::ConflictRegion,
InsertionError::Location(_) => ExplicitInsertionErrorKind::Location,
InsertionError::CavityFilling { .. } => ExplicitInsertionErrorKind::CavityFilling,
InsertionError::NeighborWiring { .. } => ExplicitInsertionErrorKind::NeighborWiring,
InsertionError::NonManifoldTopology { .. } => {
ExplicitInsertionErrorKind::NonManifoldTopology
}
InsertionError::HullExtension { .. } => ExplicitInsertionErrorKind::HullExtension,
InsertionError::DelaunayValidationFailed { .. } => {
ExplicitInsertionErrorKind::DelaunayValidationFailed
}
InsertionError::DelaunayRepairFailed { .. } => {
ExplicitInsertionErrorKind::DelaunayRepairFailed
}
InsertionError::DuplicateCoordinates { .. } => {
ExplicitInsertionErrorKind::DuplicateCoordinates
}
InsertionError::DuplicateUuid { .. } => ExplicitInsertionErrorKind::DuplicateUuid,
InsertionError::TopologyValidation(_) => ExplicitInsertionErrorKind::TopologyValidation,
InsertionError::TopologyValidationFailed { .. } => {
ExplicitInsertionErrorKind::TopologyValidationFailed
}
};
let source_kind = match &source {
InsertionError::DelaunayValidationFailed { source } => {
Some(InsertionErrorSourceKind::Delaunay(source.into()))
}
InsertionError::DelaunayRepairFailed { source, .. } => Some(
InsertionErrorSourceKind::DelaunayRepair(source.as_ref().into()),
),
InsertionError::TopologyValidation(source) => {
Some(InsertionErrorSourceKind::Tds(source.into()))
}
InsertionError::TopologyValidationFailed { source, .. } => {
Some(InsertionErrorSourceKind::Triangulation(source.into()))
}
_ => None,
};
Self {
kind,
source_kind,
message: source.to_string(),
}
}
}
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
#[non_exhaustive]
pub enum ExplicitInvariantErrorKind {
Tds,
Triangulation,
Delaunay,
}
#[must_use]
#[derive(Clone, Debug, Error, PartialEq, Eq)]
#[error("{message}")]
pub struct ExplicitInvariantError {
pub kind: ExplicitInvariantErrorKind,
pub detail: InvariantErrorSummaryDetail,
pub message: String,
}
impl From<InvariantError> for ExplicitInvariantError {
fn from(source: InvariantError) -> Self {
let kind = match &source {
InvariantError::Tds(_) => ExplicitInvariantErrorKind::Tds,
InvariantError::Triangulation(_) => ExplicitInvariantErrorKind::Triangulation,
InvariantError::Delaunay(_) => ExplicitInvariantErrorKind::Delaunay,
};
let detail = match &source {
InvariantError::Tds(source) => InvariantErrorSummaryDetail::Tds(source.into()),
InvariantError::Triangulation(source) => {
InvariantErrorSummaryDetail::Triangulation(source.into())
}
InvariantError::Delaunay(source) => {
InvariantErrorSummaryDetail::Delaunay(source.into())
}
};
Self {
kind,
detail,
message: source.to_string(),
}
}
}
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
#[non_exhaustive]
pub enum ExplicitDelaunayValidationErrorKind {
Tds,
Triangulation,
VerificationFailed,
RepairFailed,
RepairOperationFailed,
}
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
#[non_exhaustive]
pub enum ExplicitDelaunayValidationSourceKind {
Tds(TdsErrorKind),
Triangulation(TriangulationValidationErrorKind),
Repair(DelaunayRepairErrorKind),
}
#[must_use]
#[derive(Clone, Debug, Error, PartialEq, Eq)]
#[error("{message}")]
pub struct ExplicitDelaunayValidationError {
pub kind: ExplicitDelaunayValidationErrorKind,
pub source_kind: Option<ExplicitDelaunayValidationSourceKind>,
pub message: String,
}
impl From<DelaunayTriangulationValidationError> for ExplicitDelaunayValidationError {
fn from(source: DelaunayTriangulationValidationError) -> Self {
let kind = match &source {
DelaunayTriangulationValidationError::Tds(_) => {
ExplicitDelaunayValidationErrorKind::Tds
}
DelaunayTriangulationValidationError::Triangulation(_) => {
ExplicitDelaunayValidationErrorKind::Triangulation
}
DelaunayTriangulationValidationError::VerificationFailed { .. } => {
ExplicitDelaunayValidationErrorKind::VerificationFailed
}
DelaunayTriangulationValidationError::RepairFailed { .. } => {
ExplicitDelaunayValidationErrorKind::RepairFailed
}
DelaunayTriangulationValidationError::RepairOperationFailed { .. } => {
ExplicitDelaunayValidationErrorKind::RepairOperationFailed
}
};
let source_kind = match &source {
DelaunayTriangulationValidationError::Tds(source) => Some(
ExplicitDelaunayValidationSourceKind::Tds(source.as_ref().into()),
),
DelaunayTriangulationValidationError::Triangulation(source) => Some(
ExplicitDelaunayValidationSourceKind::Triangulation(source.as_ref().into()),
),
DelaunayTriangulationValidationError::RepairOperationFailed { source, .. } => Some(
ExplicitDelaunayValidationSourceKind::Repair(source.as_ref().into()),
),
DelaunayTriangulationValidationError::VerificationFailed { .. }
| DelaunayTriangulationValidationError::RepairFailed { .. } => None,
};
Self {
kind,
source_kind,
message: source.to_string(),
}
}
}
#[derive(Clone, Debug, Error, PartialEq, Eq)]
#[non_exhaustive]
pub enum ExplicitConstructionError {
#[error(
"Cell {cell_index}: vertex index {vertex_index} is out of bounds (vertex count: {bound})"
)]
IndexOutOfBounds {
cell_index: usize,
vertex_index: usize,
bound: usize,
},
#[error("Cell {cell_index}: has {actual} vertex indices, expected {expected} for a simplex")]
InvalidCellArity {
cell_index: usize,
actual: usize,
expected: usize,
},
#[error("Cell {cell_index}: contains duplicate vertex indices")]
DuplicateVertexInCell {
cell_index: usize,
},
#[error("No cells provided for explicit construction")]
EmptyCells,
#[error("Toroidal topology cannot be combined with explicit cell construction")]
IncompatibleTopology,
#[error(
"ConstructionOptions are not applicable to explicit cell construction \
and must be left at their default values"
)]
UnsupportedConstructionOptions,
#[error("Neighbor assignment failed during explicit construction: {source}")]
NeighborAssignment {
#[source]
source: ExplicitTdsError,
},
#[error("Orientation normalization failed during explicit construction: {source}")]
OrientationNormalization {
#[source]
source: ExplicitInsertionError,
},
#[error("Structural validation failed during explicit construction: {source}")]
StructuralValidation {
#[source]
source: ExplicitTdsError,
},
#[error("Topology validation failed during explicit construction: {source}")]
TopologyValidation {
#[source]
source: ExplicitInvariantError,
},
#[error("PL-manifold completion validation failed during explicit construction: {source}")]
CompletionValidation {
#[source]
source: ExplicitInvariantError,
},
#[error("Geometric nondegeneracy validation failed during explicit construction: {source}")]
GeometricNondegeneracy {
#[source]
source: ExplicitTdsError,
},
#[error("Delaunay validation failed during explicit construction: {source}")]
DelaunayValidation {
#[source]
source: ExplicitDelaunayValidationError,
},
#[error(
"Explicit non-Euclidean connectivity is not supported for {topology:?}; Level 4 quotient validation is required"
)]
UnsupportedExplicitTopology {
topology: TopologyKind,
},
}
pub struct DelaunayTriangulationBuilder<'v, T, U, const D: usize> {
vertices: &'v [Vertex<T, U, D>],
topology: Option<ToroidalSpace<D>>,
topology_guarantee: TopologyGuarantee,
construction_options: ConstructionOptions,
use_image_point_method: bool,
explicit_cells: Option<&'v [Vec<usize>]>,
global_topology: GlobalTopology<D>,
}
impl<'v, U, const D: usize> DelaunayTriangulationBuilder<'v, f64, U, D> {
#[must_use]
pub fn new(vertices: &'v [Vertex<f64, U, D>]) -> Self {
Self {
vertices,
topology: None,
topology_guarantee: TopologyGuarantee::DEFAULT,
construction_options: ConstructionOptions::default(),
use_image_point_method: false,
explicit_cells: None,
global_topology: GlobalTopology::DEFAULT,
}
}
#[must_use]
pub fn from_vertices_and_cells(
vertices: &'v [Vertex<f64, U, D>],
cells: &'v [Vec<usize>],
) -> Self {
Self::from_vertices_and_cells_generic(vertices, cells)
}
}
impl<'v, T, U, const D: usize> DelaunayTriangulationBuilder<'v, T, U, D> {
#[must_use]
pub fn from_vertices_and_cells_generic(
vertices: &'v [Vertex<T, U, D>],
cells: &'v [Vec<usize>],
) -> Self {
Self {
vertices,
topology: None,
topology_guarantee: TopologyGuarantee::DEFAULT,
construction_options: ConstructionOptions::default(),
use_image_point_method: false,
explicit_cells: Some(cells),
global_topology: GlobalTopology::DEFAULT,
}
}
#[must_use]
pub fn from_vertices(vertices: &'v [Vertex<T, U, D>]) -> Self {
Self {
vertices,
topology: None,
topology_guarantee: TopologyGuarantee::DEFAULT,
construction_options: ConstructionOptions::default(),
use_image_point_method: false,
explicit_cells: None,
global_topology: GlobalTopology::DEFAULT,
}
}
#[must_use]
pub const fn toroidal(mut self, domain: [f64; D]) -> Self {
self.topology = Some(ToroidalSpace::new(domain));
self
}
#[must_use]
pub const fn toroidal_periodic(mut self, domain: [f64; D]) -> Self {
self.topology = Some(ToroidalSpace::new(domain));
self.use_image_point_method = true;
self
}
#[must_use]
pub const fn topology_guarantee(mut self, topology_guarantee: TopologyGuarantee) -> Self {
self.topology_guarantee = topology_guarantee;
self
}
#[must_use]
pub const fn global_topology(mut self, global_topology: GlobalTopology<D>) -> Self {
self.global_topology = global_topology;
self
}
#[must_use]
pub const fn construction_options(mut self, construction_options: ConstructionOptions) -> Self {
self.construction_options = construction_options;
self
}
}
impl<T, U, const D: usize> DelaunayTriangulationBuilder<'_, T, U, D>
where
T: CoordinateScalar,
U: DataType,
{
fn validate_topology_model<M>(model: &M) -> Result<(), DelaunayTriangulationConstructionError>
where
M: GlobalTopologyModel<D>,
{
model.validate_configuration().map_err(|error| match error {
GlobalTopologyModelError::InvalidToroidalPeriod { axis, period } => {
TriangulationConstructionError::GeometricDegeneracy {
message: format!(
"Invalid toroidal domain at axis {axis}: period {period:?}; expected finite value > 0",
),
}
.into()
}
other => TriangulationConstructionError::GeometricDegeneracy {
message: format!("Invalid topology model configuration: {other}"),
}
.into(),
})
}
fn derive_periodic_facet_key(
lifted_ordered: &[(VertexKey, [i8; D])],
facet_idx: usize,
) -> Result<PeriodicFacetKey, DelaunayTriangulationConstructionError> {
periodic_facet_key_from_lifted_vertices::<D>(lifted_ordered, facet_idx).map_err(|error| {
TriangulationConstructionError::GeometricDegeneracy {
message: format!(
"Failed to derive periodic candidate facet signature for index {facet_idx}: {error}",
),
}
.into()
})
}
fn canonicalize_vertices<M>(
vertices: &[Vertex<T, U, D>],
model: &M,
) -> Result<Vec<Vertex<T, U, D>>, DelaunayTriangulationConstructionError>
where
M: GlobalTopologyModel<D>,
{
let mut out = Vec::with_capacity(vertices.len());
for (idx, v) in vertices.iter().enumerate() {
let mut canonicalized_coords = *v.point().coords();
model
.canonicalize_point_in_place(&mut canonicalized_coords)
.map_err(|error| TriangulationConstructionError::GeometricDegeneracy {
message: format!(
"Failed to canonicalize vertex {idx}: original coords {:?}; reason: {error}",
v.point().coords(),
),
})?;
let new_point = Point::new(canonicalized_coords);
let new_vertex = Vertex::new_with_uuid(new_point, v.uuid(), v.data);
out.push(new_vertex);
}
Ok(out)
}
pub fn build<V>(
self,
) -> Result<
DelaunayTriangulation<AdaptiveKernel<T>, U, V, D>,
DelaunayTriangulationConstructionError,
>
where
V: DataType,
{
self.build_with_kernel(&AdaptiveKernel::new())
}
pub fn build_with_kernel<K, V>(
self,
kernel: &K,
) -> Result<DelaunayTriangulation<K, U, V, D>, DelaunayTriangulationConstructionError>
where
K: Kernel<D, Scalar = T>,
V: DataType,
{
if let Some(cells) = self.explicit_cells {
if self.topology.is_some() {
return Err(ExplicitConstructionError::IncompatibleTopology.into());
}
if self.construction_options != ConstructionOptions::default() {
return Err(ExplicitConstructionError::UnsupportedConstructionOptions.into());
}
return Self::build_explicit(
kernel,
self.vertices,
cells,
self.topology_guarantee,
self.global_topology,
);
}
match (self.topology, self.use_image_point_method) {
(None, _) => {
let mut dt = DelaunayTriangulation::with_topology_guarantee_and_options(
kernel,
self.vertices,
self.topology_guarantee,
self.construction_options,
)?;
dt.set_global_topology(self.global_topology);
Ok(dt)
}
(Some(space), false) => {
let topology = GlobalTopology::Toroidal {
domain: space.domain,
mode: ToroidalConstructionMode::Canonicalized,
};
let topology_model = topology.model();
Self::validate_topology_model(&topology_model)?;
let canonical = Self::canonicalize_vertices(self.vertices, &topology_model)?;
let mut dt = DelaunayTriangulation::with_topology_guarantee_and_options(
kernel,
&canonical,
self.topology_guarantee,
self.construction_options,
)?;
dt.set_global_topology(topology);
Ok(dt)
}
(Some(space), true) => {
let topology = GlobalTopology::Toroidal {
domain: space.domain,
mode: ToroidalConstructionMode::PeriodicImagePoint,
};
let topology_model = topology.model();
Self::validate_topology_model(&topology_model)?;
if !topology_model.supports_periodic_facet_signatures() {
return Err(TriangulationConstructionError::GeometricDegeneracy {
message: format!(
"Topology {:?} does not support periodic facet signatures required for periodic image-point construction",
topology_model.kind(),
),
}
.into());
}
let canonical = Self::canonicalize_vertices(self.vertices, &topology_model)?;
let mut dt = Self::build_periodic(
kernel,
&canonical,
&topology_model,
self.topology_guarantee,
self.construction_options,
)?;
dt.set_global_topology(topology);
dt.tri
.normalize_and_promote_positive_orientation()
.map_err(|e| TriangulationConstructionError::GeometricDegeneracy {
message: format!(
"Failed to canonicalize periodic orientation after build: {e}",
),
})?;
dt.as_triangulation()
.validate_geometric_cell_orientation()
.map_err(|e| TriangulationConstructionError::GeometricDegeneracy {
message: format!(
"Periodic geometric orientation validation failed after build: {e}",
),
})?;
Ok(dt)
}
}
}
fn build_explicit<K, V>(
kernel: &K,
vertices: &[Vertex<T, U, D>],
cells: &[Vec<usize>],
topology_guarantee: TopologyGuarantee,
global_topology: GlobalTopology<D>,
) -> Result<DelaunayTriangulation<K, U, V, D>, DelaunayTriangulationConstructionError>
where
K: Kernel<D, Scalar = T>,
V: DataType,
{
Self::validate_explicit_cell_specs(vertices.len(), cells)?;
Self::reject_explicit_non_euclidean_topology(global_topology)?;
let vertex_count = vertices.len();
let mut tds: Tds<T, U, V, D> = Tds::empty();
let mut index_to_key = Vec::with_capacity(vertex_count);
for v in vertices {
let vk = tds
.insert_vertex_with_mapping(*v)
.map_err(TriangulationConstructionError::from)?;
index_to_key.push(vk);
}
for (cell_idx, cell_spec) in cells.iter().enumerate() {
let vertex_keys: Vec<VertexKey> =
cell_spec.iter().map(|&vi| index_to_key[vi]).collect();
let cell = Cell::new(vertex_keys, None).map_err(|e| {
TriangulationConstructionError::FailedToCreateCell {
message: format!("explicit: cell {cell_idx}: {e}"),
}
})?;
tds.insert_cell_with_mapping(cell)
.map_err(TriangulationConstructionError::from)?;
}
tds.construction_state = TriangulationConstructionState::Constructed;
tds.assign_neighbors()
.map_err(|source| ExplicitConstructionError::NeighborAssignment {
source: source.into(),
})?;
tds.assign_incident_cells().map_err(|e| {
TriangulationConstructionError::InternalInconsistency {
message: format!("explicit: incident cell assignment failed: {e}"),
}
})?;
let mut dt = DelaunayTriangulation::from_tds_with_topology_guarantee(
tds,
kernel.clone(),
topology_guarantee,
);
dt.set_global_topology(global_topology);
dt.tri
.normalize_and_promote_positive_orientation()
.map_err(
|source| ExplicitConstructionError::OrientationNormalization {
source: source.into(),
},
)?;
if let Err(e) = dt.tri.tds.validate() {
return Err(
ExplicitConstructionError::StructuralValidation { source: e.into() }.into(),
);
}
if let Err(e) = dt.tri.is_valid_topology_only() {
return Err(ExplicitConstructionError::TopologyValidation { source: e.into() }.into());
}
if let Err(e) = dt.tri.validate_at_completion() {
return Err(
ExplicitConstructionError::CompletionValidation { source: e.into() }.into(),
);
}
if let Err(e) = dt.tri.validate_geometric_nondegeneracy() {
return Err(
ExplicitConstructionError::GeometricNondegeneracy { source: e.into() }.into(),
);
}
Self::enforce_explicit_delaunay_property(&dt)?;
Ok(dt)
}
fn validate_explicit_cell_specs(
vertex_count: usize,
cells: &[Vec<usize>],
) -> Result<(), ExplicitConstructionError> {
if cells.is_empty() {
return Err(ExplicitConstructionError::EmptyCells);
}
for (cell_idx, cell_spec) in cells.iter().enumerate() {
if cell_spec.len() != D + 1 {
return Err(ExplicitConstructionError::InvalidCellArity {
cell_index: cell_idx,
actual: cell_spec.len(),
expected: D + 1,
});
}
for (i, &vi) in cell_spec.iter().enumerate() {
if vi >= vertex_count {
return Err(ExplicitConstructionError::IndexOutOfBounds {
cell_index: cell_idx,
vertex_index: vi,
bound: vertex_count,
});
}
for &vj in &cell_spec[i + 1..] {
if vi == vj {
return Err(ExplicitConstructionError::DuplicateVertexInCell {
cell_index: cell_idx,
});
}
}
}
}
Ok(())
}
fn enforce_explicit_delaunay_property<K, V>(
dt: &DelaunayTriangulation<K, U, V, D>,
) -> Result<(), DelaunayTriangulationConstructionError>
where
K: Kernel<D, Scalar = T>,
V: DataType,
{
dt.is_valid().map_err(|source| {
ExplicitConstructionError::DelaunayValidation {
source: source.into(),
}
.into()
})
}
fn reject_explicit_non_euclidean_topology(
global_topology: GlobalTopology<D>,
) -> Result<(), DelaunayTriangulationConstructionError> {
if global_topology.is_euclidean() {
return Ok(());
}
Err(ExplicitConstructionError::UnsupportedExplicitTopology {
topology: global_topology.kind(),
}
.into())
}
#[expect(
clippy::too_many_lines,
reason = "Image-point periodic DT algorithm is inherently multi-step; splitting would harm readability"
)]
fn build_periodic<K, V, M>(
kernel: &K,
canonical_vertices: &[Vertex<T, U, D>],
topology_model: &M,
topology_guarantee: TopologyGuarantee,
construction_options: ConstructionOptions,
) -> Result<DelaunayTriangulation<K, U, V, D>, DelaunayTriangulationConstructionError>
where
K: Kernel<D, Scalar = T>,
V: DataType,
M: GlobalTopologyModel<D>,
{
Self::validate_topology_model(topology_model)?;
if !topology_model.supports_periodic_facet_signatures() {
return Err(TriangulationConstructionError::GeometricDegeneracy {
message: format!(
"Topology {:?} does not support periodic facet signatures required for periodic image-point construction",
topology_model.kind(),
),
}
.into());
}
let domain = topology_model.periodic_domain().ok_or_else(|| {
TriangulationConstructionError::GeometricDegeneracy {
message: format!(
"Topology {:?} does not expose a periodic domain required for periodic image-point construction",
topology_model.kind(),
),
}
})?;
let n = canonical_vertices.len();
let min_points = 2 * D + 1;
if n < min_points {
return Err(TriangulationConstructionError::GeometricDegeneracy {
message: format!(
"Periodic {D}D triangulation requires at least {min_points} points, got {n}"
),
}
.into());
}
let three_pow_d: usize = 3_usize.pow(u32::try_from(D).expect("dimension D fits in u32"));
let zero_offset_idx = (three_pow_d - 1) / 2;
let canonical_uuids: Vec<Uuid> = canonical_vertices.iter().map(Vertex::uuid).collect();
let perturb_units = |canon_idx: usize, axis: usize| -> i64 {
let mut h = FNV_OFFSET_BASIS;
h ^= u64::try_from(canon_idx).expect("canonical index fits in u64");
h = h.wrapping_mul(FNV_PRIME);
h ^= u64::try_from(axis).expect("axis index fits in u64");
h = h.wrapping_mul(FNV_PRIME);
let span = u64::try_from(2 * MAX_OFFSET_UNITS + 1).expect("span fits in u64");
i64::try_from(h % span).expect("residue fits in i64") - MAX_OFFSET_UNITS
};
let image_jitter_units = |canon_idx: usize, axis: usize, image_idx: usize| -> i64 {
let mut h = FNV_OFFSET_BASIS;
h ^= u64::try_from(canon_idx).expect("canonical index fits in u64");
h = h.wrapping_mul(FNV_PRIME);
h ^= u64::try_from(axis).expect("axis index fits in u64");
h = h.wrapping_mul(FNV_PRIME);
h ^= u64::try_from(image_idx).expect("image index fits in u64");
h = h.wrapping_mul(FNV_PRIME);
let span = u64::try_from(2 * IMAGE_JITTER_UNITS + 1).expect("span fits in u64");
i64::try_from(h % span).expect("residue fits in i64") - IMAGE_JITTER_UNITS
};
let canonical_f64: Vec<[f64; D]> = canonical_vertices
.iter()
.enumerate()
.map(|(canon_idx, v)| {
let orig_coords = v.point().coords();
let mut coords = [0_f64; D];
for i in 0..D {
let domain_i = domain[i];
let orig = orig_coords[i]
.to_f64()
.expect("canonical coordinate is finite and convertible");
let normalized = (orig / domain_i).clamp(0.0, 1.0 - f64::EPSILON);
let u = (normalized * TWO_POW_52_F64)
.floor()
.to_i64()
.expect("grid index fits in i64");
let min_off = -u.min(MAX_OFFSET_UNITS);
let max_off = (TWO_POW_52_I64 - 1 - u).min(MAX_OFFSET_UNITS);
let off = perturb_units(canon_idx, i).clamp(min_off, max_off);
let adjusted_u = <f64 as num_traits::NumCast>::from(u + off)
.expect("adjusted grid index fits in f64");
coords[i] = (adjusted_u / TWO_POW_52_F64) * domain_i;
}
coords
})
.collect();
let mut image_uuid_to_canonical_with_offset: FastHashMap<Uuid, (Uuid, [i8; D])> =
FastHashMap::default();
let mut expanded: Vec<Vertex<T, U, D>> = Vec::with_capacity(n.saturating_mul(three_pow_d));
for k in 0..three_pow_d {
let mut offset = [0i8; D];
for (i, offset_val) in offset.iter_mut().enumerate() {
let digit =
(k / 3_usize.pow(u32::try_from(i).expect("dimension index fits in u32"))) % 3;
*offset_val = i8::try_from(digit).expect("digit is 0, 1, or 2; fits in i8") - 1;
}
let is_canonical = k == zero_offset_idx;
for (canon_idx, v) in canonical_vertices.iter().enumerate() {
let mut new_coords = [T::zero(); D];
for i in 0..D {
let shift_f64 = <f64 as From<i8>>::from(offset[i]) * domain[i];
let jitter_f64 = if is_canonical {
0.0
} else {
let jitter_units = image_jitter_units(canon_idx, i, k);
(<f64 as num_traits::NumCast>::from(jitter_units)
.expect("jitter fits in f64")
/ TWO_POW_52_F64)
* domain[i]
};
let coord_f64 = canonical_f64[canon_idx][i] + shift_f64 + jitter_f64;
new_coords[i] =
<T as num_traits::NumCast>::from(coord_f64).ok_or_else(|| {
TriangulationConstructionError::GeometricDegeneracy {
message: format!("Overflow on axis {i}: image coord {coord_f64}"),
}
})?;
}
let new_point = Point::new(new_coords);
if is_canonical {
image_uuid_to_canonical_with_offset.insert(v.uuid(), (v.uuid(), [0_i8; D]));
let canonical_v = Vertex::new_with_uuid(new_point, v.uuid(), v.data);
expanded.push(canonical_v);
} else {
let image_v: Vertex<T, U, D> = Vertex::from_point(new_point);
image_uuid_to_canonical_with_offset.insert(image_v.uuid(), (v.uuid(), offset));
expanded.push(image_v);
}
}
}
let expanded_base_options = construction_options
.with_initial_simplex_strategy(InitialSimplexStrategy::Balanced)
.without_global_repair_fallback();
let expanded_options = match construction_options.retry_policy() {
RetryPolicy::Disabled => expanded_base_options,
RetryPolicy::Shuffled { base_seed, .. }
| RetryPolicy::DebugOnlyShuffled { base_seed, .. } => expanded_base_options
.with_retry_policy(RetryPolicy::Shuffled {
attempts: NonZeroUsize::new(24).expect("literal is non-zero"),
base_seed,
}),
};
let full_dt: DelaunayTriangulation<K, U, V, D> =
match DelaunayTriangulation::with_topology_guarantee_and_options(
kernel,
&expanded,
TopologyGuarantee::Pseudomanifold,
expanded_options,
) {
Ok(dt) => dt,
Err(primary_err) if D > 2 => {
let (total_attempts, retry_seed) = match expanded_options.retry_policy() {
RetryPolicy::Disabled => (0_usize, None),
RetryPolicy::Shuffled {
attempts,
base_seed,
}
| RetryPolicy::DebugOnlyShuffled {
attempts,
base_seed,
} => (
attempts.get().saturating_mul(4).clamp(24, 256),
Some(base_seed.unwrap_or(0xA5A5_5A5A_D1E1_A1E1_u64)),
),
};
let mut built: Option<DelaunayTriangulation<K, U, V, D>> = None;
let mut last_insert_error: Option<String> = None;
let mut last_skipped_insertion: Option<String> = None;
let mut best_fallback_stats: (usize, usize, usize, usize) = (0, 0, 0, 0);
let mut insertion_order: Vec<usize> = Vec::with_capacity(expanded.len());
let canonical_start = zero_offset_idx * n;
let canonical_end = canonical_start + n;
for attempt_idx in 0..total_attempts {
insertion_order.clear();
insertion_order.extend(canonical_start..canonical_end);
insertion_order.extend(0..canonical_start);
insertion_order.extend(canonical_end..expanded.len());
if attempt_idx > 0 {
let retry_seed = retry_seed
.expect("retry_seed is only used when retry attempts are enabled");
let attempt_u64 =
u64::try_from(attempt_idx).expect("attempt index fits in u64");
let mut rng = StdRng::seed_from_u64(
retry_seed
.wrapping_add(attempt_u64.wrapping_mul(0x9E37_79B9_7F4A_7C15)),
);
let (canonical_prefix, image_suffix) = insertion_order.split_at_mut(n);
debug_assert_eq!(canonical_prefix.len(), n);
image_suffix.shuffle(&mut rng);
}
let mut candidate_dt: DelaunayTriangulation<K, U, V, D> =
DelaunayTriangulation::with_empty_kernel_and_topology_guarantee(
kernel.clone(),
TopologyGuarantee::Pseudomanifold,
);
candidate_dt.set_delaunay_repair_policy(DelaunayRepairPolicy::Never);
let mut inserted = 0_usize;
let mut skipped = 0_usize;
let mut hard_errors = 0_usize;
for (insert_idx, &source_idx) in insertion_order.iter().enumerate() {
match candidate_dt.insert_with_statistics(expanded[source_idx]) {
Ok((InsertionOutcome::Inserted { .. }, _stats)) => {
inserted = inserted.saturating_add(1);
}
Ok((InsertionOutcome::Skipped { error }, _stats)) => {
skipped = skipped.saturating_add(1);
last_skipped_insertion = Some(format!(
"attempt={attempt_idx} insert_idx={insert_idx} source_idx={source_idx}: {error}",
));
}
Err(err) => {
hard_errors = hard_errors.saturating_add(1);
last_insert_error = Some(format!(
"attempt={attempt_idx} insert_idx={insert_idx} source_idx={source_idx}: {err}",
));
}
}
}
let canonical_present = canonical_uuids
.iter()
.filter(|uuid| candidate_dt.tds().vertex_key_from_uuid(uuid).is_some())
.count();
if canonical_present > best_fallback_stats.0
|| (canonical_present == best_fallback_stats.0
&& inserted > best_fallback_stats.1)
{
best_fallback_stats =
(canonical_present, inserted, skipped, hard_errors);
}
if canonical_present == n
&& candidate_dt.number_of_cells() > 0
&& candidate_dt.tds().is_valid().is_ok()
{
built = Some(candidate_dt);
break;
}
}
if let Some(dt) = built {
dt
} else {
let canonical_vertex_uuid_sample: Vec<Uuid> = canonical_vertices
.iter()
.take(3)
.map(Vertex::uuid)
.collect();
return Err(TriangulationConstructionError::GeometricDegeneracy {
message: format!(
"Periodic expanded DT construction failed (no fallback): canonical_vertices_len={}, canonical_vertex_uuid_sample={canonical_vertex_uuid_sample:?}, primary_err={primary_err}, last_insert_error={:?}, last_skipped_insertion={:?}, best_fallback_stats(canonical_present,inserted,skipped,hard_errors)={:?}, topology_guarantee={topology_guarantee:?}, construction_options={construction_options:?}",
canonical_vertices.len(),
last_insert_error,
last_skipped_insertion,
best_fallback_stats,
),
}
.into());
}
}
Err(err) => return Err(err),
};
let tds_ref = full_dt.tds();
let Some(central_keys) = canonical_uuids
.iter()
.map(|uuid| tds_ref.vertex_key_from_uuid(uuid))
.collect::<Option<Vec<_>>>()
else {
return Err(TriangulationConstructionError::GeometricDegeneracy {
message: format!(
"Periodic expanded DT is missing at least one canonical vertex: canonical_vertices_len={}",
canonical_uuids.len()
),
}
.into());
};
let central_key_set: VertexKeySet = central_keys.into_iter().collect();
let mut vertex_key_to_lifted: FastHashMap<VertexKey, (VertexKey, [i8; D])> =
FastHashMap::default();
for vk in tds_ref.vertex_keys() {
let Some(vertex) = tds_ref.vertex(vk) else {
continue;
};
let Some((canonical_uuid, offset)) =
image_uuid_to_canonical_with_offset.get(&vertex.uuid())
else {
continue;
};
let Some(canonical_key) = tds_ref.vertex_key_from_uuid(canonical_uuid) else {
continue;
};
vertex_key_to_lifted.insert(vk, (canonical_key, *offset));
}
let normalize_cell_lifted = |cell_key: CellKey| -> Option<Vec<(VertexKey, [i8; D])>> {
let cell = tds_ref.cell(cell_key)?;
let mut lifted: Vec<(VertexKey, [i8; D])> = cell
.vertices()
.iter()
.map(|vk| vertex_key_to_lifted.get(vk).copied())
.collect::<Option<Vec<_>>>()?;
let mut canonical_keys: Vec<VertexKey> = lifted.iter().map(|(ck, _)| *ck).collect();
canonical_keys.sort_unstable();
canonical_keys.dedup();
if canonical_keys.len() != D + 1 {
return None;
}
let (anchor_idx, _) = lifted.iter().enumerate().min_by_key(|(_, (ck, _))| *ck)?;
let anchor_offset = lifted[anchor_idx].1;
for (_, offset) in &mut lifted {
for axis in 0..D {
offset[axis] -= anchor_offset[axis];
}
}
Some(lifted)
};
let cell_barycenter_in_fundamental_domain = |cell_key: CellKey| -> Option<bool> {
let cell = tds_ref.cell(cell_key)?;
let mut sums = [0.0_f64; D];
for vk in cell.vertices() {
let vertex = tds_ref.vertex(*vk)?;
let coords = vertex.point().coords();
for (axis, sum) in sums.iter_mut().enumerate() {
*sum += coords[axis].to_f64()?;
}
}
let denom = <f64 as num_traits::NumCast>::from(D + 1)
.expect("simplex vertex count fits in f64 for D");
for (axis, sum) in sums.iter().enumerate() {
let bary = *sum / denom;
let period = domain[axis];
if !(bary >= 0.0 && bary < period) {
return Some(false);
}
}
Some(true)
};
let mut candidates_by_symbolic: FastHashMap<SymbolicSignature<D>, PeriodicCandidate<D>> =
FastHashMap::default();
for ck in tds_ref.cell_keys() {
let Some(lifted_vertices) = normalize_cell_lifted(ck) else {
continue;
};
let in_domain = cell_barycenter_in_fundamental_domain(ck).unwrap_or(false);
let mut symbolic_signature = lifted_vertices.clone();
symbolic_signature.sort_unstable();
let lifted_ordered = lifted_vertices.clone();
let mut periodic_facets: Vec<PeriodicFacetKey> = Vec::with_capacity(D + 1);
for facet_idx in 0..=D {
periodic_facets.push(Self::derive_periodic_facet_key(&lifted_ordered, facet_idx)?);
}
if let Some(existing) = candidates_by_symbolic.get_mut(&symbolic_signature) {
if in_domain {
existing.3 = true;
}
} else {
candidates_by_symbolic.insert(
symbolic_signature.clone(),
(
symbolic_signature,
lifted_ordered,
periodic_facets,
in_domain,
),
);
}
}
let mut candidates: Vec<PeriodicCandidate<D>> =
candidates_by_symbolic.into_values().collect();
if candidates.is_empty() {
return Err(TriangulationConstructionError::GeometricDegeneracy {
message: "No quotient periodic cells found in full image DT".to_owned(),
}
.into());
}
candidates.sort_by(|a, b| b.3.cmp(&a.3).then_with(|| a.0.cmp(&b.0)));
let (search_attempts, search_seed) = match construction_options.retry_policy() {
RetryPolicy::Disabled => (1_usize, 0xD1CE_0B5E_2100_0001_u64),
RetryPolicy::Shuffled {
attempts,
base_seed,
}
| RetryPolicy::DebugOnlyShuffled {
attempts,
base_seed,
} => (
attempts
.get()
.saturating_add(1)
.saturating_mul(512)
.clamp(512, 4096),
base_seed.unwrap_or(0xD1CE_0B5E_2100_0001_u64),
),
};
let mut best_selected: Vec<bool> = Vec::new();
let mut best_boundary_count = usize::MAX;
let mut best_selected_count = 0_usize;
let mut best_coverage_count = 0_usize;
let mut best_abs_chi = i64::MAX;
if D == 2 {
let target_faces = central_key_set.len().saturating_mul(2);
let mut edge_to_index: FastHashMap<PeriodicFacetKey, usize> = FastHashMap::default();
let mut candidate_edges: Vec<[usize; 3]> = Vec::with_capacity(candidates.len());
let mut candidate_in_domain: Vec<bool> = Vec::with_capacity(candidates.len());
for candidate in &candidates {
let mut edge_indices = [0usize; 3];
for (slot, edge_key) in candidate.2.iter().enumerate() {
let next_index = edge_to_index.len();
let edge_index = *edge_to_index.entry(*edge_key).or_insert(next_index);
edge_indices[slot] = edge_index;
}
candidate_edges.push(edge_indices);
candidate_in_domain.push(candidate.3);
}
let exact_search_node_limit = candidate_edges
.len()
.saturating_mul(edge_to_index.len().max(1))
.saturating_mul(512)
.clamp(100_000, 5_000_000);
if let Some(exact_selected) = search_closed_2d_selection(
&candidate_edges,
&candidate_in_domain,
target_faces,
edge_to_index.len(),
exact_search_node_limit,
) {
best_selected_count = exact_selected
.iter()
.filter(|&&is_selected| is_selected)
.count();
best_coverage_count = central_key_set.len();
best_boundary_count = 0;
best_abs_chi = 0;
best_selected = exact_selected;
}
}
if best_selected.is_empty() {
let base_order: Vec<usize> = (0..candidates.len()).collect();
for attempt_idx in 0..search_attempts {
let mut order = base_order.clone();
if attempt_idx > 0 {
let attempt_u64 =
u64::try_from(attempt_idx).expect("attempt index fits in u64");
let mut rng = StdRng::seed_from_u64(
search_seed.wrapping_add(attempt_u64.wrapping_mul(0x9E37_79B9_7F4A_7C15)),
);
order.shuffle(&mut rng);
}
order.sort_by(|a, b| candidates[*b].3.cmp(&candidates[*a].3));
let mut selected = vec![false; candidates.len()];
let mut facet_counts: FastHashMap<PeriodicFacetKey, u8> = FastHashMap::default();
for idx in order.iter().copied() {
let candidate_facets = &candidates[idx].2;
if candidate_facets
.iter()
.any(|facet| facet_counts.get(facet).copied().unwrap_or(0) >= 2)
{
continue;
}
selected[idx] = true;
for facet in candidate_facets {
*facet_counts.entry(*facet).or_insert(0) += 1;
}
}
let mut improved = true;
while improved {
improved = false;
for idx in order.iter().copied() {
if selected[idx] {
continue;
}
let candidate_facets = &candidates[idx].2;
if candidate_facets
.iter()
.any(|facet| facet_counts.get(facet).copied().unwrap_or(0) >= 2)
{
continue;
}
let boundary_delta: i32 = candidate_facets
.iter()
.map(
|facet| match facet_counts.get(facet).copied().unwrap_or(0) {
0 => 1,
1 => -1,
_ => 0,
},
)
.sum();
if boundary_delta < 0 {
selected[idx] = true;
for facet in candidate_facets {
*facet_counts.entry(*facet).or_insert(0) += 1;
}
improved = true;
}
}
}
loop {
let mut best_move: Option<(bool, usize, i32)> = None;
for idx in order.iter().copied() {
let candidate_facets = &candidates[idx].2;
if selected[idx] {
let boundary_delta: i32 = candidate_facets
.iter()
.map(
|facet| match facet_counts.get(facet).copied().unwrap_or(0) {
1 => -1,
2 => 1,
_ => 0,
},
)
.sum();
if boundary_delta < 0
&& best_move
.is_none_or(|(_, _, best_delta)| boundary_delta < best_delta)
{
best_move = Some((false, idx, boundary_delta));
}
} else {
if candidate_facets
.iter()
.any(|facet| facet_counts.get(facet).copied().unwrap_or(0) >= 2)
{
continue;
}
let boundary_delta: i32 = candidate_facets
.iter()
.map(
|facet| match facet_counts.get(facet).copied().unwrap_or(0) {
0 => 1,
1 => -1,
_ => 0,
},
)
.sum();
if boundary_delta < 0
&& best_move
.is_none_or(|(_, _, best_delta)| boundary_delta < best_delta)
{
best_move = Some((true, idx, boundary_delta));
}
}
}
let Some((is_add, idx, _)) = best_move else {
break;
};
let candidate_facets = &candidates[idx].2;
if is_add {
selected[idx] = true;
for facet in candidate_facets {
*facet_counts.entry(*facet).or_insert(0) += 1;
}
} else {
selected[idx] = false;
for facet in candidate_facets {
if let Some(count) = facet_counts.get_mut(facet) {
*count -= 1;
if *count == 0 {
facet_counts.remove(facet);
}
}
}
}
}
let boundary_count = facet_counts.values().filter(|&&count| count == 1).count();
let selected_count = selected.iter().filter(|&&is_selected| is_selected).count();
let mut covered: VertexKeySet = VertexKeySet::default();
for (idx, is_selected) in selected.iter().copied().enumerate() {
if !is_selected {
continue;
}
for (vertex_key, _) in &candidates[idx].1 {
covered.insert(*vertex_key);
}
}
let coverage_count = covered.len();
let abs_chi = if D == 2 {
let v_count =
i64::try_from(central_key_set.len()).expect("vertex count fits in i64");
let e_count =
i64::try_from(facet_counts.len()).expect("edge/facet count fits in i64");
let f_count = i64::try_from(selected_count).expect("cell count fits in i64");
(v_count - e_count + f_count).abs()
} else {
0
};
if boundary_count < best_boundary_count
|| (boundary_count == best_boundary_count
&& (if D == 2 {
abs_chi < best_abs_chi
|| (abs_chi == best_abs_chi && selected_count > best_selected_count)
} else {
coverage_count > best_coverage_count
|| (coverage_count == best_coverage_count
&& selected_count > best_selected_count)
}))
{
best_boundary_count = boundary_count;
best_selected_count = selected_count;
best_coverage_count = coverage_count;
best_abs_chi = abs_chi;
best_selected = selected;
}
let best_has_full_canonical_coverage = best_coverage_count == central_key_set.len();
if best_boundary_count == 0
&& ((D == 2 && best_abs_chi == 0)
|| (D > 2 && best_has_full_canonical_coverage))
{
break;
}
}
}
if best_selected.is_empty() {
return Err(TriangulationConstructionError::GeometricDegeneracy {
message: "Periodic quotient selection failed to pick any candidate cells"
.to_owned(),
}
.into());
}
if D == 2 && best_boundary_count > 0 {
return Err(TriangulationConstructionError::GeometricDegeneracy {
message: format!(
"Periodic quotient selection left {best_boundary_count} boundary facets after {search_attempts} attempts (full_vertices={}, full_cells={}, canonical_vertices={}, candidates={}, selected_cells={})",
tds_ref.number_of_vertices(),
tds_ref.number_of_cells(),
central_key_set.len(),
candidates.len(),
best_selected_count,
),
}
.into());
}
if D == 2 && best_abs_chi != 0 {
return Err(TriangulationConstructionError::GeometricDegeneracy {
message: format!(
"Periodic quotient selection could not reach χ=0 in 2D (best |χ|={best_abs_chi}) after {search_attempts} attempts",
),
}
.into());
}
let has_full_canonical_coverage = best_coverage_count == central_key_set.len();
if D > 2 && !has_full_canonical_coverage {
return Err(TriangulationConstructionError::GeometricDegeneracy {
message: format!(
"Periodic quotient selection covered only {} of {} canonical vertices in {D}D",
best_coverage_count,
central_key_set.len(),
),
}
.into());
}
if D > 2 {
let mut selected_by_canonical: FastHashMap<Vec<VertexKey>, usize> =
FastHashMap::default();
let mut dedup_selected = vec![false; candidates.len()];
for (idx, is_selected) in best_selected.iter().copied().enumerate() {
if !is_selected {
continue;
}
let mut canonical_keys: Vec<VertexKey> =
candidates[idx].1.iter().map(|(vk, _)| *vk).collect();
canonical_keys.sort_unstable();
if let Some(existing_idx) = selected_by_canonical.get(&canonical_keys).copied() {
let existing_in_domain = candidates[existing_idx].3;
let candidate_in_domain = candidates[idx].3;
let should_replace = (!existing_in_domain && candidate_in_domain)
|| (existing_in_domain == candidate_in_domain
&& candidates[idx].0 < candidates[existing_idx].0);
if should_replace {
dedup_selected[existing_idx] = false;
dedup_selected[idx] = true;
selected_by_canonical.insert(canonical_keys, idx);
}
} else {
dedup_selected[idx] = true;
selected_by_canonical.insert(canonical_keys, idx);
}
}
best_selected = dedup_selected;
}
let mut representative_lifted_by_symbolic: FastHashMap<
SymbolicSignature<D>,
SymbolicSignature<D>,
> = FastHashMap::default();
for (idx, is_selected) in best_selected.iter().copied().enumerate() {
if !is_selected {
continue;
}
let (symbolic_signature, lifted_ordered, _, _) = &candidates[idx];
representative_lifted_by_symbolic
.insert(symbolic_signature.clone(), lifted_ordered.clone());
}
let mut tds_mut = tds_ref.clone();
let all_cells: Vec<CellKey> = tds_mut.cell_keys().collect();
tds_mut.remove_cells_by_keys(&all_cells);
let image_vertex_keys: Vec<VertexKey> = tds_mut
.vertex_keys()
.filter(|vk| !central_key_set.contains(vk))
.collect();
for &vk in &image_vertex_keys {
tds_mut.remove_vertex(vk).map_err(|e| {
TriangulationConstructionError::InternalInconsistency {
message: format!("Failed to remove image vertex: {e}"),
}
})?;
}
let mut signatures_sorted: Vec<Vec<(VertexKey, [i8; D])>> =
representative_lifted_by_symbolic.keys().cloned().collect();
signatures_sorted.sort_unstable();
let mut inserted_cell_keys: Vec<CellKey> = Vec::with_capacity(signatures_sorted.len());
let mut rep_lifted_by_key: FastHashMap<CellKey, Vec<(VertexKey, [i8; D])>> =
FastHashMap::default();
for signature in signatures_sorted {
let Some(lifted_vertices) = representative_lifted_by_symbolic.get(&signature) else {
continue;
};
let canonical_vertex_keys: Vec<VertexKey> =
lifted_vertices.iter().map(|(ck, _)| *ck).collect();
let mut cell = Cell::new(canonical_vertex_keys, None).map_err(|e| {
TriangulationConstructionError::GeometricDegeneracy {
message: format!("Failed to create quotient periodic cell: {e}"),
}
})?;
let offsets: PeriodicOffsetBuffer<D> =
lifted_vertices.iter().map(|(_, offset)| *offset).collect();
cell.set_periodic_vertex_offsets(offsets).map_err(|e| {
TriangulationConstructionError::GeometricDegeneracy {
message: format!("Failed to set quotient periodic offsets: {e}"),
}
})?;
let ck = tds_mut
.insert_cell_with_mapping_trusted_vertices(cell)
.map_err(|e| TriangulationConstructionError::GeometricDegeneracy {
message: format!("Failed to insert quotient periodic cell: {e}"),
})?;
inserted_cell_keys.push(ck);
rep_lifted_by_key.insert(ck, lifted_vertices.clone());
}
if inserted_cell_keys.is_empty() {
return Err(TriangulationConstructionError::GeometricDegeneracy {
message: "No cells survived periodic quotient reconstruction".to_owned(),
}
.into());
}
let mut periodic_facet_counts: FastHashMap<PeriodicFacetKey, usize> =
FastHashMap::default();
for lifted in rep_lifted_by_key.values() {
for facet_idx in 0..=D {
let periodic_facet_key = Self::derive_periodic_facet_key(lifted, facet_idx)?;
*periodic_facet_counts.entry(periodic_facet_key).or_insert(0) += 1;
}
}
let overloaded_facets: Vec<(PeriodicFacetKey, usize)> = periodic_facet_counts
.into_iter()
.filter(|(_, count)| *count > 2)
.collect();
if !overloaded_facets.is_empty() {
return Err(TriangulationConstructionError::GeometricDegeneracy {
message: format!(
"Periodic quotient selection overcounts periodic facets ({} overloaded); selected_cells={}, sample={:?}",
overloaded_facets.len(),
rep_lifted_by_key.len(),
overloaded_facets.iter().take(4).collect::<Vec<_>>(),
),
}
.into());
}
let mut neighbor_updates: FastHashMap<CellKey, Vec<Option<CellKey>>> = inserted_cell_keys
.iter()
.copied()
.map(|ck| (ck, vec![None; D + 1]))
.collect();
let mut facet_occurrences: FastHashMap<PeriodicFacetKey, Vec<(CellKey, usize)>> =
FastHashMap::default();
for &rep_ck in &inserted_cell_keys {
let Some(lifted) = rep_lifted_by_key.get(&rep_ck) else {
continue;
};
for facet_idx in 0..=D {
let sig = Self::derive_periodic_facet_key(lifted, facet_idx)?;
facet_occurrences
.entry(sig)
.or_default()
.push((rep_ck, facet_idx));
}
}
for (_facet_sig, occurrences) in facet_occurrences {
match occurrences.as_slice() {
[(a_ck, a_idx), (b_ck, b_idx)] => {
let a_lifted = rep_lifted_by_key.get(a_ck).ok_or_else(|| {
TriangulationConstructionError::InternalInconsistency {
message: format!(
"missing lifted representative for quotient cell {a_ck:?}"
),
}
})?;
let b_lifted = rep_lifted_by_key.get(b_ck).ok_or_else(|| {
TriangulationConstructionError::InternalInconsistency {
message: format!(
"missing lifted representative for quotient cell {b_ck:?}"
),
}
})?;
let shares_all_canonical_vertices = a_lifted
.iter()
.zip(b_lifted.iter())
.all(|((a_vk, _), (b_vk, _))| a_vk == b_vk);
if shares_all_canonical_vertices {
return Err(TriangulationConstructionError::GeometricDegeneracy {
message: format!(
"Periodic quotient produced distinct cells with identical canonical vertices across a shared facet: {a_ck:?}[{a_idx}] <-> {b_ck:?}[{b_idx}]",
),
}
.into());
}
neighbor_updates.get_mut(a_ck).ok_or_else(|| {
TriangulationConstructionError::InternalInconsistency {
message: format!("missing neighbor vector for quotient cell {a_ck:?}"),
}
})?[*a_idx] = Some(*b_ck);
neighbor_updates.get_mut(b_ck).ok_or_else(|| {
TriangulationConstructionError::InternalInconsistency {
message: format!("missing neighbor vector for quotient cell {b_ck:?}"),
}
})?[*b_idx] = Some(*a_ck);
}
[(a_ck, a_idx)] => {
neighbor_updates.get_mut(a_ck).ok_or_else(|| {
TriangulationConstructionError::InternalInconsistency {
message: format!("missing neighbor vector for quotient cell {a_ck:?}"),
}
})?[*a_idx] = Some(*a_ck);
}
_ => {
return Err(TriangulationConstructionError::GeometricDegeneracy {
message: format!(
"Periodic quotient facet signature has {} occurrences (expected 1 or 2): {occurrences:?}",
occurrences.len()
),
}
.into());
}
}
}
let unmatched_count = neighbor_updates
.values()
.flat_map(|n| n.iter())
.filter(|n| n.is_none())
.count();
if unmatched_count > 0 {
return Err(TriangulationConstructionError::GeometricDegeneracy {
message: format!(
"Periodic quotient reconstruction left {unmatched_count} unmatched neighbor slots"
),
}
.into());
}
for &ck in &inserted_cell_keys {
let neighbors = neighbor_updates.remove(&ck).ok_or_else(|| {
TriangulationConstructionError::InternalInconsistency {
message: format!("missing neighbor vector for inserted quotient cell {ck:?}"),
}
})?;
tds_mut.set_neighbors_by_key(ck, &neighbors).map_err(|e| {
TriangulationConstructionError::InternalInconsistency {
message: format!("set_neighbors_by_key failed for {ck:?}: {e}"),
}
})?;
}
if let Err(_error) = tds_mut.normalize_coherent_orientation() {
#[cfg(debug_assertions)]
tracing::debug!(
?_error,
"periodic quotient: skipping coherent-orientation normalization failure"
);
}
tds_mut.assign_incident_cells().map_err(|e| {
TriangulationConstructionError::InternalInconsistency {
message: format!("assign_incident_cells failed: {e}"),
}
})?;
if let Err(e) = tds_mut.is_valid() {
return Err(TriangulationConstructionError::GeometricDegeneracy {
message: format!("Periodic quotient TDS invalid before return: {e}"),
}
.into());
}
Ok(DelaunayTriangulation::from_tds_with_topology_guarantee(
tds_mut,
kernel.clone(),
topology_guarantee,
))
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::core::algorithms::flips::DelaunayRepairError;
use crate::core::algorithms::incremental_insertion::{
CavityFillingError, DelaunayRepairFailureContext, HullExtensionReason, NeighborWiringError,
};
use crate::core::algorithms::locate::{ConflictError, LocateError};
use crate::core::cell::CellValidationError;
use crate::core::facet::FacetError;
use crate::core::tds::{
DelaunayValidationErrorKind, EntityKind, GeometricError, NeighborValidationError,
};
use crate::core::triangulation::TriangulationValidationError;
use crate::core::util::uuid::UuidValidationError;
use crate::core::vertex::VertexBuilder;
use crate::core::vertex::VertexValidationError;
use crate::geometry::kernel::RobustKernel;
use crate::topology::traits::global_topology_model::{
EuclideanModel, GlobalTopologyModel, GlobalTopologyModelError, ToroidalModel,
};
use crate::topology::traits::topological_space::{
GlobalTopology, TopologyKind, ToroidalConstructionMode,
};
use crate::triangulation::delaunay::{
DelaunayConstructionFailure, DelaunayRepairOperation, InsertionOrderStrategy,
};
use crate::vertex;
use approx::assert_relative_eq;
use slotmap::{Key, KeyData};
#[derive(Clone, Copy, Debug)]
struct ValidationFailureModel;
impl GlobalTopologyModel<2> for ValidationFailureModel {
fn kind(&self) -> TopologyKind {
TopologyKind::Euclidean
}
fn allows_boundary(&self) -> bool {
true
}
fn validate_configuration(&self) -> Result<(), GlobalTopologyModelError> {
Err(GlobalTopologyModelError::NonFiniteCoordinate {
axis: 0,
value: f64::NAN,
})
}
fn canonicalize_point_in_place<T>(
&self,
_coords: &mut [T; 2],
) -> Result<(), GlobalTopologyModelError>
where
T: CoordinateScalar,
{
Ok(())
}
fn lift_for_orientation<T>(
&self,
coords: [T; 2],
periodic_offset: Option<[i8; 2]>,
) -> Result<[T; 2], GlobalTopologyModelError>
where
T: CoordinateScalar,
{
if periodic_offset.is_some() {
return Err(GlobalTopologyModelError::PeriodicOffsetsUnsupported {
kind: TopologyKind::Euclidean,
});
}
Ok(coords)
}
}
#[test]
fn explicit_error_summaries_preserve_nested_source_kinds() {
let tds = ExplicitTdsError::from(TdsError::DuplicateCells {
message: "same vertex set appears twice".to_string(),
});
assert_eq!(tds.kind, ExplicitTdsErrorKind::DuplicateCells);
assert!(tds.message.contains("Duplicate cells"));
let insertion = ExplicitInsertionError::from(InsertionError::TopologyValidation(
TdsError::InconsistentDataStructure {
message: "dangling neighbor".to_string(),
},
));
assert_eq!(
insertion.source_kind,
Some(InsertionErrorSourceKind::Tds(
TdsErrorKind::InconsistentDataStructure,
))
);
let invariant = ExplicitInvariantError::from(InvariantError::Triangulation(
TriangulationValidationError::Disconnected { cell_count: 2 },
));
assert_eq!(invariant.kind, ExplicitInvariantErrorKind::Triangulation);
assert_eq!(
invariant.detail,
InvariantErrorSummaryDetail::Triangulation(
TriangulationValidationErrorKind::Disconnected,
)
);
let delaunay = ExplicitDelaunayValidationError::from(
DelaunayTriangulationValidationError::RepairOperationFailed {
operation: DelaunayRepairOperation::VertexRemoval,
source: Box::new(DelaunayRepairError::HeuristicRebuildFailed {
message: "rebuild failed".to_string(),
}),
},
);
assert_eq!(
delaunay.source_kind,
Some(ExplicitDelaunayValidationSourceKind::Repair(
DelaunayRepairErrorKind::HeuristicRebuildFailed,
))
);
let delaunay_tds = ExplicitDelaunayValidationError::from(
DelaunayTriangulationValidationError::from(TdsError::InconsistentDataStructure {
message: "dangling cell".to_string(),
}),
);
assert_eq!(delaunay_tds.kind, ExplicitDelaunayValidationErrorKind::Tds);
assert_eq!(
delaunay_tds.source_kind,
Some(ExplicitDelaunayValidationSourceKind::Tds(
TdsErrorKind::InconsistentDataStructure,
))
);
let delaunay_topology =
ExplicitDelaunayValidationError::from(DelaunayTriangulationValidationError::from(
TriangulationValidationError::Disconnected { cell_count: 2 },
));
assert_eq!(
delaunay_topology.kind,
ExplicitDelaunayValidationErrorKind::Triangulation,
);
assert_eq!(
delaunay_topology.source_kind,
Some(ExplicitDelaunayValidationSourceKind::Triangulation(
TriangulationValidationErrorKind::Disconnected,
))
);
}
fn assert_explicit_tds_error_kind(source: TdsError, expected_kind: ExplicitTdsErrorKind) {
let summary = ExplicitTdsError::from(source);
assert_eq!(summary.kind, expected_kind);
assert!(!summary.message.is_empty());
}
#[test]
fn explicit_tds_error_preserves_validation_error_kinds() {
let cell_key = CellKey::from(KeyData::from_ffi(1));
let other_cell_key = CellKey::from(KeyData::from_ffi(2));
let vertex_key = VertexKey::from(KeyData::from_ffi(3));
let uuid = Uuid::new_v4();
assert_explicit_tds_error_kind(
TdsError::InvalidVertex {
vertex_id: uuid,
source: VertexValidationError::InvalidUuid {
source: UuidValidationError::NilUuid,
},
},
ExplicitTdsErrorKind::InvalidVertex,
);
assert_explicit_tds_error_kind(
TdsError::InvalidCell {
cell_id: uuid,
source: CellValidationError::DuplicateVertices,
},
ExplicitTdsErrorKind::InvalidCell,
);
assert_explicit_tds_error_kind(
TdsError::InvalidNeighbors {
reason: NeighborValidationError::Other {
message: "neighbor invariant failed".to_string(),
},
},
ExplicitTdsErrorKind::InvalidNeighbors,
);
assert_explicit_tds_error_kind(
TdsError::OrientationViolation {
cell1_key: cell_key,
cell1_uuid: uuid,
cell2_key: other_cell_key,
cell2_uuid: Uuid::new_v4(),
cell1_facet_index: 0,
cell2_facet_index: 1,
facet_vertices: vec![vertex_key],
cell2_facet_vertices: vec![vertex_key],
observed_odd_permutation: false,
expected_odd_permutation: true,
},
ExplicitTdsErrorKind::OrientationViolation,
);
assert_explicit_tds_error_kind(
TdsError::Geometric(GeometricError::DegenerateOrientation {
message: "zero determinant".to_string(),
}),
ExplicitTdsErrorKind::Geometric,
);
assert_explicit_tds_error_kind(
TdsError::FacetError(FacetError::InvalidFacetIndex {
index: 4,
facet_count: 4,
}),
ExplicitTdsErrorKind::FacetError,
);
assert_explicit_tds_error_kind(
TdsError::DuplicateCoordinatesInCell {
cell_id: uuid,
message: "two vertices share coordinates".to_string(),
},
ExplicitTdsErrorKind::DuplicateCoordinatesInCell,
);
}
#[test]
fn explicit_tds_error_preserves_lookup_and_operation_error_kinds() {
let cell_key = CellKey::from(KeyData::from_ffi(1));
let vertex_key = VertexKey::from(KeyData::from_ffi(3));
let uuid = Uuid::new_v4();
assert_explicit_tds_error_kind(
TdsError::DuplicateCells {
message: "duplicate cell vertex set".to_string(),
},
ExplicitTdsErrorKind::DuplicateCells,
);
assert_explicit_tds_error_kind(
TdsError::FacetSharingViolation {
facet_key: 42,
existing_incident_count: 2,
attempted_incident_count: 3,
max_incident_count: 2,
candidate_cell_uuid: uuid,
candidate_facet_index: 1,
},
ExplicitTdsErrorKind::FacetSharingViolation,
);
assert_explicit_tds_error_kind(
TdsError::FailedToCreateCell {
message: "cell validation failed".to_string(),
},
ExplicitTdsErrorKind::FailedToCreateCell,
);
assert_explicit_tds_error_kind(
TdsError::NotNeighbors {
cell1: uuid,
cell2: Uuid::new_v4(),
},
ExplicitTdsErrorKind::NotNeighbors,
);
assert_explicit_tds_error_kind(
TdsError::MappingInconsistency {
entity: EntityKind::Cell,
message: "uuid mapping was stale".to_string(),
},
ExplicitTdsErrorKind::MappingInconsistency,
);
assert_explicit_tds_error_kind(
TdsError::VertexKeyRetrievalFailed {
cell_id: uuid,
message: "cell vertices unavailable".to_string(),
},
ExplicitTdsErrorKind::VertexKeyRetrievalFailed,
);
assert_explicit_tds_error_kind(
TdsError::CellNotFound {
cell_key,
context: "cell lookup".to_string(),
},
ExplicitTdsErrorKind::CellNotFound,
);
assert_explicit_tds_error_kind(
TdsError::VertexNotFound {
vertex_key,
context: "vertex lookup".to_string(),
},
ExplicitTdsErrorKind::VertexNotFound,
);
assert_explicit_tds_error_kind(
TdsError::DimensionMismatch {
expected: 4,
actual: 3,
context: "simplex arity".to_string(),
},
ExplicitTdsErrorKind::DimensionMismatch,
);
assert_explicit_tds_error_kind(
TdsError::IndexOutOfBounds {
index: 4,
bound: 4,
context: "facet index".to_string(),
},
ExplicitTdsErrorKind::IndexOutOfBounds,
);
assert_explicit_tds_error_kind(
TdsError::InconsistentDataStructure {
message: "dangling neighbor".to_string(),
},
ExplicitTdsErrorKind::InconsistentDataStructure,
);
}
fn assert_explicit_insertion_error(
source: InsertionError,
expected_kind: ExplicitInsertionErrorKind,
expected_source_kind: Option<InsertionErrorSourceKind>,
) {
let summary = ExplicitInsertionError::from(source);
assert_eq!(summary.kind, expected_kind);
assert_eq!(summary.source_kind, expected_source_kind);
assert!(!summary.message.is_empty());
}
#[test]
fn explicit_insertion_error_preserves_stage_kinds_without_nested_sources() {
let cell_key = CellKey::from(KeyData::from_ffi(1));
let uuid = Uuid::new_v4();
assert_explicit_insertion_error(
InsertionError::ConflictRegion(ConflictError::InvalidStartCell { cell_key }),
ExplicitInsertionErrorKind::ConflictRegion,
None,
);
assert_explicit_insertion_error(
InsertionError::Location(LocateError::EmptyTriangulation),
ExplicitInsertionErrorKind::Location,
None,
);
assert_explicit_insertion_error(
InsertionError::CavityFilling {
reason: CavityFillingError::MissingBoundaryCell { cell_key },
},
ExplicitInsertionErrorKind::CavityFilling,
None,
);
assert_explicit_insertion_error(
InsertionError::NeighborWiring {
reason: NeighborWiringError::MissingCell { cell_key },
},
ExplicitInsertionErrorKind::NeighborWiring,
None,
);
assert_explicit_insertion_error(
InsertionError::NonManifoldTopology {
facet_hash: 0xabc,
cell_count: 3,
},
ExplicitInsertionErrorKind::NonManifoldTopology,
None,
);
assert_explicit_insertion_error(
InsertionError::HullExtension {
reason: HullExtensionReason::NoVisibleFacets,
},
ExplicitInsertionErrorKind::HullExtension,
None,
);
assert_explicit_insertion_error(
InsertionError::DuplicateCoordinates {
coordinates: "[0.0, 0.0]".to_string(),
},
ExplicitInsertionErrorKind::DuplicateCoordinates,
None,
);
assert_explicit_insertion_error(
InsertionError::DuplicateUuid {
entity: EntityKind::Vertex,
uuid,
},
ExplicitInsertionErrorKind::DuplicateUuid,
None,
);
}
#[test]
fn explicit_insertion_error_preserves_nested_validation_source_kinds() {
assert_explicit_insertion_error(
InsertionError::DelaunayValidationFailed {
source: DelaunayTriangulationValidationError::VerificationFailed {
message: "non-Delaunay facet".to_string(),
},
},
ExplicitInsertionErrorKind::DelaunayValidationFailed,
Some(InsertionErrorSourceKind::Delaunay(
DelaunayValidationErrorKind::VerificationFailed,
)),
);
assert_explicit_insertion_error(
InsertionError::DelaunayRepairFailed {
source: Box::new(DelaunayRepairError::PostconditionFailed {
message: "remaining violation".to_string(),
}),
context: DelaunayRepairFailureContext::LocalRepair,
},
ExplicitInsertionErrorKind::DelaunayRepairFailed,
Some(InsertionErrorSourceKind::DelaunayRepair(
DelaunayRepairErrorKind::PostconditionFailed,
)),
);
assert_explicit_insertion_error(
InsertionError::TopologyValidation(TdsError::Geometric(
GeometricError::DegenerateOrientation {
message: "zero determinant".to_string(),
},
)),
ExplicitInsertionErrorKind::TopologyValidation,
Some(InsertionErrorSourceKind::Tds(TdsErrorKind::Geometric)),
);
assert_explicit_insertion_error(
InsertionError::TopologyValidationFailed {
source: TriangulationValidationError::RidgeLinkNotManifold {
ridge_key: 0xdef,
link_vertex_count: 4,
link_edge_count: 2,
max_degree: 3,
degree_one_vertices: 1,
connected: false,
},
message: "scoped topology validation".to_string(),
},
ExplicitInsertionErrorKind::TopologyValidationFailed,
Some(InsertionErrorSourceKind::Triangulation(
TriangulationValidationErrorKind::RidgeLinkNotManifold,
)),
);
}
#[derive(Clone, Copy, Debug)]
struct CanonicalizationFailureModel;
impl GlobalTopologyModel<2> for CanonicalizationFailureModel {
fn kind(&self) -> TopologyKind {
TopologyKind::Euclidean
}
fn allows_boundary(&self) -> bool {
true
}
fn validate_configuration(&self) -> Result<(), GlobalTopologyModelError> {
Ok(())
}
fn canonicalize_point_in_place<T>(
&self,
_coords: &mut [T; 2],
) -> Result<(), GlobalTopologyModelError>
where
T: CoordinateScalar,
{
Err(GlobalTopologyModelError::NonFiniteCoordinate {
axis: 0,
value: f64::NAN,
})
}
fn lift_for_orientation<T>(
&self,
coords: [T; 2],
periodic_offset: Option<[i8; 2]>,
) -> Result<[T; 2], GlobalTopologyModelError>
where
T: CoordinateScalar,
{
if periodic_offset.is_some() {
return Err(GlobalTopologyModelError::PeriodicOffsetsUnsupported {
kind: TopologyKind::Euclidean,
});
}
Ok(coords)
}
}
#[derive(Clone, Copy, Debug)]
struct MissingPeriodicDomainModel;
impl GlobalTopologyModel<2> for MissingPeriodicDomainModel {
fn kind(&self) -> TopologyKind {
TopologyKind::Toroidal
}
fn allows_boundary(&self) -> bool {
false
}
fn validate_configuration(&self) -> Result<(), GlobalTopologyModelError> {
Ok(())
}
fn canonicalize_point_in_place<T>(
&self,
_coords: &mut [T; 2],
) -> Result<(), GlobalTopologyModelError>
where
T: CoordinateScalar,
{
Ok(())
}
fn lift_for_orientation<T>(
&self,
coords: [T; 2],
_periodic_offset: Option<[i8; 2]>,
) -> Result<[T; 2], GlobalTopologyModelError>
where
T: CoordinateScalar,
{
Ok(coords)
}
fn supports_periodic_facet_signatures(&self) -> bool {
true
}
fn periodic_domain(&self) -> Option<&[f64; 2]> {
None
}
}
#[test]
fn test_builder_euclidean_2d() {
let vertices = vec![
vertex!([0.0, 0.0]),
vertex!([1.0, 0.0]),
vertex!([0.0, 1.0]),
];
let dt = DelaunayTriangulationBuilder::new(&vertices)
.build::<()>()
.unwrap();
assert_eq!(dt.number_of_vertices(), 3);
assert_eq!(dt.dim(), 2);
assert!(dt.as_triangulation().validate().is_ok());
}
#[test]
fn test_builder_euclidean_3d() {
let vertices = vec![
vertex!([0.0, 0.0, 0.0]),
vertex!([1.0, 0.0, 0.0]),
vertex!([0.0, 1.0, 0.0]),
vertex!([0.0, 0.0, 1.0]),
];
let dt = DelaunayTriangulationBuilder::new(&vertices)
.build::<()>()
.unwrap();
assert_eq!(dt.number_of_vertices(), 4);
assert!(dt.validate().is_ok());
}
#[test]
fn test_builder_topology_guarantee_propagated() {
let vertices = vec![
vertex!([0.0, 0.0]),
vertex!([1.0, 0.0]),
vertex!([0.0, 1.0]),
];
let dt = DelaunayTriangulationBuilder::new(&vertices)
.topology_guarantee(TopologyGuarantee::Pseudomanifold)
.build::<()>()
.unwrap();
assert_eq!(dt.topology_guarantee(), TopologyGuarantee::Pseudomanifold);
}
#[test]
fn test_builder_custom_options_propagated() {
let vertices = vec![
vertex!([0.0, 0.0]),
vertex!([1.0, 0.0]),
vertex!([0.0, 1.0]),
];
let opts =
ConstructionOptions::default().with_insertion_order(InsertionOrderStrategy::Input);
let dt = DelaunayTriangulationBuilder::new(&vertices)
.construction_options(opts)
.build::<()>()
.unwrap();
assert_eq!(dt.number_of_vertices(), 3);
}
#[test]
fn test_builder_toroidal_canonicalizes_out_of_domain_vertices() {
let vertices = vec![
vertex!([0.2, 0.3]), vertex!([1.8, 0.1]), vertex!([0.5, 0.7]), vertex!([-0.4, 0.9]), ];
let dt = DelaunayTriangulationBuilder::new(&vertices)
.toroidal([1.0, 1.0])
.build::<()>()
.unwrap();
for (_, v) in dt.vertices() {
let c = v.point().coords();
assert!(c[0] >= 0.0 && c[0] < 1.0, "x = {} not in [0, 1)", c[0]);
assert!(c[1] >= 0.0 && c[1] < 1.0, "y = {} not in [0, 1)", c[1]);
}
assert_eq!(dt.number_of_vertices(), 4);
}
#[test]
fn test_builder_toroidal_in_domain_vertices_unchanged() {
let vertices = vec![
vertex!([0.1, 0.2]),
vertex!([0.8, 0.3]),
vertex!([0.4, 0.9]),
];
let dt = DelaunayTriangulationBuilder::new(&vertices)
.toroidal([1.0, 1.0])
.build::<()>()
.unwrap();
for (_, v) in dt.vertices() {
let c = v.point().coords();
assert!(c[0] >= 0.0 && c[0] < 1.0);
assert!(c[1] >= 0.0 && c[1] < 1.0);
}
}
#[test]
fn test_builder_toroidal_build_succeeds_2d() {
let vertices = vec![
vertex!([0.2, 0.3]),
vertex!([0.8, 0.1]),
vertex!([0.5, 0.7]),
vertex!([0.1, 0.9]),
];
let dt = DelaunayTriangulationBuilder::new(&vertices)
.toroidal([1.0, 1.0])
.build::<()>()
.unwrap();
assert_eq!(dt.number_of_vertices(), 4);
assert_eq!(dt.dim(), 2);
assert!(dt.as_triangulation().validate().is_ok());
assert!(matches!(
dt.global_topology(),
GlobalTopology::Toroidal {
mode: ToroidalConstructionMode::Canonicalized,
..
}
));
}
#[test]
fn test_builder_toroidal_build_out_of_domain_input_2d() {
let vertices = vec![
vertex!([2.2, 3.3]), vertex!([-0.2, 1.1]), vertex!([1.5, 0.7]), vertex!([-0.9, 2.9]), ];
let dt = DelaunayTriangulationBuilder::new(&vertices)
.toroidal([1.0, 1.0])
.build::<()>()
.unwrap();
assert_eq!(dt.number_of_vertices(), 4);
assert_eq!(dt.dim(), 2);
assert!(dt.as_triangulation().validate().is_ok());
}
#[test]
fn test_builder_toroidal_non_finite_coordinate_is_error() {
let vertices = vec![
VertexBuilder::<f64, (), 2>::default()
.point(Point::new([0.2_f64, 0.3]))
.build()
.unwrap(),
VertexBuilder::<f64, (), 2>::default()
.point(Point::new([f64::NAN, 0.1]))
.build()
.unwrap(),
VertexBuilder::<f64, (), 2>::default()
.point(Point::new([0.5_f64, 0.7]))
.build()
.unwrap(),
];
let result = DelaunayTriangulationBuilder::new(&vertices)
.toroidal([1.0, 1.0])
.build::<()>();
assert!(result.is_err());
}
#[test]
fn test_builder_toroidal_invalid_domain_is_error() {
let vertices = vec![
vertex!([0.2, 0.3]),
vertex!([0.8, 0.1]),
vertex!([0.5, 0.7]),
];
let result = DelaunayTriangulationBuilder::new(&vertices)
.toroidal([0.0, 1.0])
.build::<()>();
let err = result.expect_err("zero period should be rejected");
assert!(format!("{err}").contains("Invalid toroidal domain"));
}
#[test]
fn test_builder_toroidal_periodic_invalid_domain_is_error() {
let vertices = vec![
vertex!([0.1, 0.2]),
vertex!([0.4, 0.7]),
vertex!([0.7, 0.3]),
vertex!([0.2, 0.9]),
vertex!([0.8, 0.6]),
vertex!([0.5, 0.1]),
vertex!([0.3, 0.5]),
];
let result = DelaunayTriangulationBuilder::new(&vertices)
.toroidal_periodic([1.0, 0.0])
.build::<()>();
let err = result.expect_err("zero period should be rejected");
assert!(format!("{err}").contains("Invalid toroidal domain"));
}
#[test]
fn test_builder_toroidal_periodic_2d_smoke() {
let vertices = vec![
vertex!([0.1_f64, 0.2]),
vertex!([0.4, 0.7]),
vertex!([0.7, 0.3]),
vertex!([0.2, 0.9]),
vertex!([0.8, 0.6]),
vertex!([0.5, 0.1]),
vertex!([0.3, 0.5]),
];
let n = vertices.len();
let kernel = RobustKernel::new();
let dt = DelaunayTriangulationBuilder::new(&vertices)
.toroidal_periodic([1.0, 1.0])
.build_with_kernel::<_, ()>(&kernel)
.unwrap();
assert_eq!(dt.number_of_vertices(), n);
assert!(dt.tds().is_valid().is_ok());
assert!(matches!(
dt.global_topology(),
GlobalTopology::Toroidal {
mode: ToroidalConstructionMode::PeriodicImagePoint,
..
}
));
}
#[test]
fn test_builder_toroidal_idempotent_on_canonical_input() {
let vertices = vec![
vertex!([0.1, 0.2]),
vertex!([0.8, 0.3]),
vertex!([0.4, 0.9]),
];
let dt_euclidean = DelaunayTriangulationBuilder::new(&vertices)
.build::<()>()
.unwrap();
let dt_toroidal = DelaunayTriangulationBuilder::new(&vertices)
.toroidal([1.0, 1.0])
.build::<()>()
.unwrap();
assert_eq!(
dt_euclidean.number_of_vertices(),
dt_toroidal.number_of_vertices()
);
assert_eq!(
dt_euclidean.number_of_cells(),
dt_toroidal.number_of_cells()
);
}
#[test]
fn test_builder_from_vertices_preserves_vertex_data() {
let vertices: Vec<Vertex<f64, i32, 2>> = vec![
VertexBuilder::default()
.point(Point::new([0.2_f64, 0.3]))
.data(1_i32)
.build()
.unwrap(),
VertexBuilder::default()
.point(Point::new([1.8_f64, 0.1])) .data(2_i32)
.build()
.unwrap(),
VertexBuilder::default()
.point(Point::new([0.5_f64, 0.7]))
.data(3_i32)
.build()
.unwrap(),
];
let dt = DelaunayTriangulationBuilder::new(&vertices)
.toroidal([1.0, 1.0])
.build::<()>()
.unwrap();
assert_eq!(dt.number_of_vertices(), 3);
for (_, v) in dt.vertices() {
let c = v.point().coords();
assert!(c[0] >= 0.0 && c[0] < 1.0);
assert!(c[1] >= 0.0 && c[1] < 1.0);
}
let mut data: Vec<i32> = dt.vertices().filter_map(|(_, v)| v.data).collect();
data.sort_unstable();
assert_eq!(data, vec![1, 2, 3]);
}
#[test]
fn test_builder_with_robust_kernel() {
let vertices = vec![
vertex!([0.0, 0.0, 0.0]),
vertex!([1.0, 0.0, 0.0]),
vertex!([0.0, 1.0, 0.0]),
vertex!([0.0, 0.0, 1.0]),
];
let kernel = RobustKernel::<f64>::new();
let dt = DelaunayTriangulationBuilder::new(&vertices)
.build_with_kernel::<_, ()>(&kernel)
.unwrap();
assert_eq!(dt.number_of_vertices(), 4);
assert!(dt.validate().is_ok());
}
#[test]
fn test_validate_topology_model_accepts_valid_toroidal() {
let model = ToroidalModel::<2>::new([2.0, 3.0], ToroidalConstructionMode::Canonicalized);
let result = DelaunayTriangulationBuilder::<f64, (), 2>::validate_topology_model(&model);
assert!(result.is_ok());
}
#[test]
fn test_derive_periodic_facet_key_happy_path_matches_core_derivation() {
let lifted_ordered = vec![
(VertexKey::null(), [0_i8, 0_i8]),
(VertexKey::null(), [1_i8, 0_i8]),
(VertexKey::null(), [0_i8, 1_i8]),
];
let expected = periodic_facet_key_from_lifted_vertices::<2>(&lifted_ordered, 1).unwrap();
let derived = DelaunayTriangulationBuilder::<f64, (), 2>::derive_periodic_facet_key(
&lifted_ordered,
1,
)
.unwrap();
assert_eq!(derived, expected);
}
#[test]
fn test_derive_periodic_facet_key_maps_errors_to_geometric_degeneracy() {
let lifted_ordered = vec![
(VertexKey::null(), [-128_i8, 0_i8]),
(VertexKey::null(), [0_i8, 0_i8]),
(VertexKey::null(), [127_i8, 0_i8]),
];
let err = DelaunayTriangulationBuilder::<f64, (), 2>::derive_periodic_facet_key(
&lifted_ordered,
1,
)
.unwrap_err();
match err {
DelaunayTriangulationConstructionError::Triangulation(
DelaunayConstructionFailure::GeometricDegeneracy { message },
) => {
assert!(
message.contains(
"Failed to derive periodic candidate facet signature for index 1"
),
"unexpected message prefix: {message}"
);
assert!(
message.contains("out of encodable range"),
"expected wrapped derivation detail in message: {message}"
);
}
other => panic!("expected GeometricDegeneracy mapping, got: {other:?}"),
}
}
#[test]
fn test_validate_topology_model_rejects_zero_period() {
let model = ToroidalModel::<2>::new([0.0, 3.0], ToroidalConstructionMode::Canonicalized);
let result = DelaunayTriangulationBuilder::<f64, (), 2>::validate_topology_model(&model);
assert!(result.is_err());
let err_str = format!("{}", result.unwrap_err());
assert!(
err_str.contains("Invalid toroidal domain"),
"Error message should mention invalid toroidal domain: {err_str}"
);
assert!(
err_str.contains("axis 0"),
"Error message should mention axis: {err_str}"
);
}
#[test]
fn test_validate_topology_model_rejects_negative_period() {
let model =
ToroidalModel::<3>::new([2.0, -1.0, 3.0], ToroidalConstructionMode::Canonicalized);
let result = DelaunayTriangulationBuilder::<f64, (), 3>::validate_topology_model(&model);
assert!(result.is_err());
let err_str = format!("{}", result.unwrap_err());
assert!(err_str.contains("Invalid toroidal domain"));
assert!(err_str.contains("axis 1"));
}
#[test]
fn test_validate_topology_model_rejects_infinite_period() {
let model = ToroidalModel::<2>::new(
[f64::INFINITY, 3.0],
ToroidalConstructionMode::Canonicalized,
);
let result = DelaunayTriangulationBuilder::<f64, (), 2>::validate_topology_model(&model);
assert!(result.is_err());
let err_str = format!("{}", result.unwrap_err());
assert!(err_str.contains("Invalid toroidal domain"));
}
#[test]
fn test_validate_topology_model_rejects_nan_period() {
let model =
ToroidalModel::<2>::new([f64::NAN, 3.0], ToroidalConstructionMode::Canonicalized);
let result = DelaunayTriangulationBuilder::<f64, (), 2>::validate_topology_model(&model);
assert!(result.is_err());
let err_str = format!("{}", result.unwrap_err());
assert!(err_str.contains("Invalid toroidal domain"));
}
#[test]
fn test_validate_topology_model_accepts_euclidean() {
let model = EuclideanModel;
let result = DelaunayTriangulationBuilder::<f64, (), 2>::validate_topology_model(&model);
assert!(result.is_ok());
}
#[test]
fn test_validate_topology_model_maps_non_period_errors() {
let result = DelaunayTriangulationBuilder::<f64, (), 2>::validate_topology_model(
&ValidationFailureModel,
);
let err = result.expect_err("non-period validation failure should be mapped");
let err_str = err.to_string();
assert!(err_str.contains("Invalid topology model configuration"));
}
#[test]
fn test_canonicalize_vertices_preserves_uuids() {
let vertices = vec![
vertex!([2.5, 3.7]),
vertex!([1.8, -0.5]),
vertex!([0.5, 0.7]),
];
let original_uuids: Vec<_> = vertices.iter().map(Vertex::uuid).collect();
let model = ToroidalModel::<2>::new([2.0, 3.0], ToroidalConstructionMode::Canonicalized);
let canonical =
DelaunayTriangulationBuilder::<f64, (), 2>::canonicalize_vertices(&vertices, &model)
.unwrap();
assert_eq!(canonical.len(), vertices.len());
let canonical_uuids: Vec<_> = canonical.iter().map(Vertex::uuid).collect();
assert_eq!(canonical_uuids, original_uuids);
}
#[test]
fn test_canonicalize_vertices_preserves_data() {
let vertices: Vec<Vertex<f64, i32, 2>> = vec![
VertexBuilder::default()
.point(Point::new([2.5_f64, 3.7]))
.data(10_i32)
.build()
.unwrap(),
VertexBuilder::default()
.point(Point::new([1.8_f64, -0.5]))
.data(20_i32)
.build()
.unwrap(),
VertexBuilder::default()
.point(Point::new([0.5_f64, 0.7]))
.data(30_i32)
.build()
.unwrap(),
];
let model = ToroidalModel::<2>::new([2.0, 3.0], ToroidalConstructionMode::Canonicalized);
let canonical =
DelaunayTriangulationBuilder::<f64, i32, 2>::canonicalize_vertices(&vertices, &model)
.unwrap();
assert_eq!(canonical.len(), vertices.len());
for (orig, canon) in vertices.iter().zip(canonical.iter()) {
assert_eq!(orig.data, canon.data);
}
}
#[test]
fn test_canonicalize_vertices_transforms_coordinates() {
let vertices = vec![
vertex!([2.5, 3.7]), vertex!([1.8, -0.5]), vertex!([0.3, 0.2]), ];
let model = ToroidalModel::<2>::new([2.0, 3.0], ToroidalConstructionMode::Canonicalized);
let canonical =
DelaunayTriangulationBuilder::<f64, (), 2>::canonicalize_vertices(&vertices, &model)
.unwrap();
assert_eq!(canonical.len(), 3);
assert_relative_eq!(canonical[0].point().coords()[0], 0.5);
assert_relative_eq!(canonical[0].point().coords()[1], 0.7);
assert_relative_eq!(canonical[1].point().coords()[0], 1.8);
assert_relative_eq!(canonical[1].point().coords()[1], 2.5);
assert_relative_eq!(canonical[2].point().coords()[0], 0.3);
assert_relative_eq!(canonical[2].point().coords()[1], 0.2);
}
#[test]
fn test_canonicalize_vertices_includes_vertex_context_on_error() {
let vertices = vec![vertex!([0.25_f64, 0.75_f64]), vertex!([0.9_f64, 0.1_f64])];
let result = DelaunayTriangulationBuilder::<f64, (), 2>::canonicalize_vertices(
&vertices,
&CanonicalizationFailureModel,
);
let err = result.expect_err("canonicalization failure should be reported");
let err_str = err.to_string();
assert!(err_str.contains("Failed to canonicalize vertex 0"));
assert!(err_str.contains("reason"));
}
#[test]
fn test_build_periodic_requires_periodic_domain() {
let kernel = AdaptiveKernel::new();
let canonical_vertices = vec![
vertex!([0.1_f64, 0.1_f64]),
vertex!([0.9_f64, 0.2_f64]),
vertex!([0.2_f64, 0.8_f64]),
vertex!([0.7_f64, 0.9_f64]),
vertex!([0.5_f64, 0.4_f64]),
];
let result = DelaunayTriangulationBuilder::<f64, (), 2>::build_periodic::<_, (), _>(
&kernel,
&canonical_vertices,
&MissingPeriodicDomainModel,
TopologyGuarantee::default(),
ConstructionOptions::default(),
);
let err = result.expect_err("missing periodic domain must fail");
assert!(err.to_string().contains(
"does not expose a periodic domain required for periodic image-point construction"
));
}
#[test]
fn test_canonicalize_vertices_euclidean_identity() {
let vertices = vec![
vertex!([1.5, 2.5]),
vertex!([3.7, 4.2]),
vertex!([-1.0, -2.0]),
];
let model = EuclideanModel;
let canonical =
DelaunayTriangulationBuilder::<f64, (), 2>::canonicalize_vertices(&vertices, &model)
.unwrap();
assert_eq!(canonical.len(), vertices.len());
for (orig, canon) in vertices.iter().zip(canonical.iter()) {
assert_relative_eq!(orig.point().coords()[0], canon.point().coords()[0]);
assert_relative_eq!(orig.point().coords()[1], canon.point().coords()[1]);
}
}
#[test]
fn test_canonicalize_vertices_propagates_nan_error() {
let vertices = vec![
VertexBuilder::default()
.point(Point::new([0.5_f64, 0.5]))
.build()
.unwrap(),
VertexBuilder::default()
.point(Point::new([f64::NAN, 0.5]))
.build()
.unwrap(),
VertexBuilder::default()
.point(Point::new([0.3_f64, 0.2]))
.build()
.unwrap(),
];
let model = ToroidalModel::<2>::new([2.0, 3.0], ToroidalConstructionMode::Canonicalized);
let result =
DelaunayTriangulationBuilder::<f64, (), 2>::canonicalize_vertices(&vertices, &model);
assert!(result.is_err());
let err_str = format!("{}", result.unwrap_err());
assert!(
err_str.contains("Failed to canonicalize vertex"),
"Error should mention canonicalization failure: {err_str}"
);
assert!(
err_str.contains("vertex 1"),
"Error should mention vertex index: {err_str}"
);
}
#[test]
fn test_canonicalize_vertices_propagates_infinity_error() {
let vertices = vec![
VertexBuilder::default()
.point(Point::new([0.5_f64, 0.5]))
.build()
.unwrap(),
VertexBuilder::default()
.point(Point::new([0.3_f64, 0.2]))
.build()
.unwrap(),
VertexBuilder::default()
.point(Point::new([f64::INFINITY, 0.5]))
.build()
.unwrap(),
];
let model = ToroidalModel::<2>::new([2.0, 3.0], ToroidalConstructionMode::Canonicalized);
let result =
DelaunayTriangulationBuilder::<f64, (), 2>::canonicalize_vertices(&vertices, &model);
assert!(result.is_err());
let err_str = format!("{}", result.unwrap_err());
assert!(err_str.contains("Failed to canonicalize vertex"));
assert!(err_str.contains("vertex 2"));
}
#[test]
fn test_canonicalize_vertices_includes_original_coords_in_error() {
let vertices = vec![
VertexBuilder::default()
.point(Point::new([f64::NAN, 1.5_f64]))
.build()
.unwrap(),
];
let model = ToroidalModel::<2>::new([2.0, 3.0], ToroidalConstructionMode::Canonicalized);
let result =
DelaunayTriangulationBuilder::<f64, (), 2>::canonicalize_vertices(&vertices, &model);
assert!(result.is_err());
let err_str = format!("{}", result.unwrap_err());
assert!(
err_str.contains("original coords"),
"Error should mention original coords: {err_str}"
);
}
}