bherbruckandClaude Opus 5.5
296683229d
* WIP deep rename (not for push) * WIP2 * WIP3 * WIP4 wxs fix * rustfmt * Rename the remaining identifiers in merged crash-report code Co-Authored-By: Claude Opus 5.5 <[email protected]> Claude-Session: https://claude.ai/code/session_01RjGw4gA5ZvKrU8MQP6sBgT --------- Co-authored-by: Claude Opus 5.5 <[email protected]>
904 lines
39 KiB
Rust
904 lines
39 KiB
Rust
//! Independent current-geometry intervals for uniform prisms.
|
|
//! These are not mesh measurements or certificates of construction history.
|
|
use crate::body::from_p3;
|
|
use solidcraft_geom::{Plane, Seg2, Vec2, Vec3};
|
|
use truck_modeling as mt;
|
|
|
|
#[derive(Clone, Copy, Debug)]
|
|
struct Interval {
|
|
lo: f64,
|
|
hi: f64,
|
|
}
|
|
|
|
impl Interval {
|
|
fn enclosing(value: f64, error: f64) -> Option<Self> {
|
|
if !value.is_finite() || !error.is_finite() || error < 0. {
|
|
return None;
|
|
}
|
|
// The explicit extra arithmetic allowance also covers subtraction and
|
|
// accumulation; do not round an uncertainty inward at its endpoints.
|
|
let padding = 16. * f64::EPSILON * (value.abs() + error + 1.);
|
|
Some(Self { lo: value - error - padding, hi: value + error + padding })
|
|
}
|
|
|
|
fn add(self, other: Self) -> Option<Self> {
|
|
Self::enclosing((self.lo + self.hi + other.lo + other.hi) * 0.5, (self.hi - self.lo + other.hi - other.lo) * 0.5)
|
|
}
|
|
}
|
|
|
|
/// Signed Green integral of a fitted segment with the *whole-curve* spatial
|
|
/// uncertainty. Angular/longitudinal monotonicity is mandatory: a curve in a
|
|
/// thin corridor could otherwise retrace it arbitrarily often.
|
|
fn green_interval(segment: &Seg2, tolerance: f64) -> Option<Interval> {
|
|
let (length, a, b) = (segment.length(), segment.start(), segment.end());
|
|
let extra = match *segment {
|
|
Seg2::Line { .. } => 0.,
|
|
Seg2::Arc { center, radius, .. } => center.len() + radius,
|
|
_ => return None,
|
|
};
|
|
let scale = length + a.len() + b.len() + extra;
|
|
// A monotone line differs from its chord by a strip; a radial-monotone
|
|
// arc contributes 1/2 integral r(theta)^2 dtheta. The endpoint wedges,
|
|
// coordinate uncertainty, and angular endpoint uncertainty are included
|
|
// by this conservative bound (rather than treating the fitted arc exact).
|
|
let error = 8. * tolerance * scale + 64. * tolerance * tolerance + 1024. * f64::EPSILON * (1. + scale * scale);
|
|
Interval::enclosing(segment.area_term(), error)
|
|
}
|
|
|
|
struct CapLoop {
|
|
segments: Vec<Seg2>,
|
|
area: Interval,
|
|
}
|
|
|
|
/// Cross-product cancellation belongs to the surface representation, not
|
|
/// just the later incidence dot product. Reject ill-conditioned plane bases.
|
|
fn plane_normal(p: &mt::Plane) -> Option<(Vec3, Vec3)> {
|
|
let u = p.u_axis();
|
|
let v = p.v_axis();
|
|
let u = Vec3::new(u.x, u.y, u.z);
|
|
let v = Vec3::new(v.x, v.y, v.z);
|
|
let cross = u.cross(v);
|
|
let length = cross.len();
|
|
if !length.is_finite() || length <= 1e-6 * u.len() * v.len() {
|
|
return None;
|
|
}
|
|
let n = cross.normalized()?;
|
|
let ce = Vec3::new((u.y * v.z).abs() + (u.z * v.y).abs(), (u.z * v.x).abs() + (u.x * v.z).abs(), (u.x * v.y).abs() + (u.y * v.x).abs())
|
|
* (256. * f64::EPSILON);
|
|
let ne = ce * (1. / length) + Vec3::new(n.x.abs(), n.y.abs(), n.z.abs()) * (ce.len() / length);
|
|
(ne.x.is_finite() && ne.y.is_finite() && ne.z.is_finite()).then_some((n, ne))
|
|
}
|
|
|
|
/// A whole positive control-net bound for the wall's projected variation.
|
|
/// A sampled normal or a vertex-only test is deliberately insufficient.
|
|
fn wall_error(surface: &mt::Surface, d: Vec3, height: f64, work: &mut usize) -> Option<f64> {
|
|
if let mt::Surface::Plane(p) = surface {
|
|
let (n, ne) = plane_normal(p)?;
|
|
let terms = (n.x * d.x).abs() + (n.y * d.y).abs() + (n.z * d.z).abs();
|
|
let error = ne.x * d.x.abs() + ne.y * d.y.abs() + ne.z * d.z.abs() + 1024. * f64::EPSILON * terms;
|
|
return Some((n.dot(d).abs() + error) * height * 2.);
|
|
}
|
|
if let mt::Surface::RevolutedCurve(p) = surface {
|
|
let r = p.entity();
|
|
let mt::Curve::Line(line) = r.entity_curve() else { return None };
|
|
let axis = r.axis();
|
|
let meridian = line.1 - line.0;
|
|
let a = Vec3::new(axis.x, axis.y, axis.z).normalized()?;
|
|
let m = Vec3::new(meridian.x, meridian.y, meridian.z).normalized()?;
|
|
let placed = *p.transform() * mt::Vector4::new(axis.x, axis.y, axis.z, 0.);
|
|
if placed.w != 0. {
|
|
return None;
|
|
}
|
|
let pa = Vec3::new(placed.x, placed.y, placed.z).normalized()?;
|
|
return Some((a.cross(m).len() + pa.cross(d).len()) * height * 4.);
|
|
}
|
|
let (rows, cols) = match surface {
|
|
mt::Surface::BSplineSurface(s) => (s.control_points().len(), s.control_points().first()?.len()),
|
|
mt::Surface::NurbsSurface(s) => (s.control_points().len(), s.control_points().first()?.len()),
|
|
_ => return None,
|
|
};
|
|
let count = rows.checked_mul(cols)?;
|
|
let rectangular = match surface {
|
|
mt::Surface::BSplineSurface(s) => s.control_points().iter().all(|r| r.len() == cols),
|
|
mt::Surface::NurbsSurface(s) => s.control_points().iter().all(|r| r.len() == cols),
|
|
_ => false,
|
|
};
|
|
if rows < 2 || cols < 2 || count > 16_384 || !rectangular {
|
|
return None;
|
|
}
|
|
*work = work.checked_sub(count.checked_mul(4)?)?;
|
|
let points: Vec<Vec<(Vec3, f64)>> = match surface {
|
|
mt::Surface::BSplineSurface(s) => s.control_points().iter().map(|r| r.iter().map(|p| (from_p3(*p), 1.)).collect()).collect(),
|
|
mt::Surface::NurbsSurface(s) => {
|
|
s.control_points().iter().map(|r| r.iter().map(|p| (Vec3::new(p.x / p.w, p.y / p.w, p.z / p.w), p.w)).collect()).collect()
|
|
}
|
|
_ => return None,
|
|
};
|
|
if points.iter().any(|r| r.len() != cols) || points.iter().flatten().any(|(p, w)| !p.is_finite() || !w.is_finite() || *w <= 0.) {
|
|
return None;
|
|
}
|
|
let max_w = points.iter().flatten().map(|p| p.1).fold(0., f64::max);
|
|
let min_w = points.iter().flatten().map(|p| p.1 / max_w).fold(1., f64::min);
|
|
if min_w < 1e-6 {
|
|
return None;
|
|
}
|
|
let magnitude = points.iter().flatten().map(|p| p.0.len()).fold(1., f64::max);
|
|
let mut best = f64::INFINITY;
|
|
for transpose in [false, true] {
|
|
let at = |i: usize, j: usize| if transpose { points.get(j)?.get(i).copied() } else { points.get(i)?.get(j).copied() };
|
|
let (groups, length) = if transpose { (cols, rows) } else { (rows, cols) };
|
|
let mut variation: f64 = 0.;
|
|
let mut rank_error: f64 = 0.;
|
|
let mut monotone = true;
|
|
let mut axial_sign = 0.;
|
|
for i in 0..groups {
|
|
let (base, bw) = at(i, 0)?;
|
|
let mut previous = base.dot(d);
|
|
for j in 0..length {
|
|
let (point, w) = at(i, j)?;
|
|
let step = point.dot(d) - previous;
|
|
previous = point.dot(d);
|
|
if step.abs() > 2048. * f64::EPSILON * magnitude {
|
|
if axial_sign * step < 0. {
|
|
monotone = false;
|
|
}
|
|
axial_sign = step.signum();
|
|
}
|
|
variation = variation.max((point - base).cross(d).len());
|
|
let lhs = w / max_w * (at(0, 0)?.1 / max_w);
|
|
let rhs = bw / max_w * (at(0, j)?.1 / max_w);
|
|
rank_error = rank_error.max((lhs - rhs).abs());
|
|
}
|
|
if (at(i, length - 1)?.0 - base).dot(d).abs() < height * 0.5 {
|
|
monotone = false;
|
|
}
|
|
}
|
|
if !monotone || axial_sign == 0. {
|
|
continue;
|
|
}
|
|
// Positive weights turn the control hull into a positional bound.
|
|
// Rank-one weight error perturbs the normalized convex coefficients;
|
|
// the positive minimum weight prevents an unbounded denominator.
|
|
let error = 2. * variation + 8. * magnitude * rank_error / min_w.powi(2) + 2048. * f64::EPSILON * magnitude / min_w.powi(2);
|
|
best = best.min(error);
|
|
}
|
|
best.is_finite().then_some(best)
|
|
}
|
|
|
|
fn segment_same(a: &Seg2, b: &Seg2, tolerance: f64) -> bool {
|
|
match (*a, *b) {
|
|
(Seg2::Line { a, b }, Seg2::Line { a: c, b: d }) => {
|
|
(a.dist(c) <= tolerance && b.dist(d) <= tolerance) || (a.dist(d) <= tolerance && b.dist(c) <= tolerance)
|
|
}
|
|
(Seg2::Arc { center, radius, start, sweep }, Seg2::Arc { center: c, radius: r, start: t, sweep: s }) => {
|
|
center.dist(c) <= tolerance
|
|
&& (radius - r).abs() <= tolerance
|
|
&& (sweep.abs() - s.abs()).abs() * radius <= tolerance
|
|
&& ((Vec2::new(start.cos(), start.sin()) * radius).dist(Vec2::new(t.cos(), t.sin()) * r) <= tolerance
|
|
|| (Vec2::new(start.cos(), start.sin()) * radius).dist(Vec2::new((t + s).cos(), (t + s).sin()) * r) <= tolerance)
|
|
}
|
|
_ => false,
|
|
}
|
|
}
|
|
|
|
fn wall_supports_cap(surface: &mt::Surface, edge: &mt::Edge, plane: &Plane, levels: (f64, f64), tol: f64, work: &mut usize) -> Option<bool> {
|
|
use mt::{BoundedCurve, ParametricCurve};
|
|
let cap = crate::prism_curve::segment(edge, plane, tol, work)?;
|
|
if let mt::Surface::Plane(p) = surface {
|
|
let (n, ne) = plane_normal(p)?;
|
|
let origin = from_p3(p.origin());
|
|
for s in cap {
|
|
let Seg2::Line { a, b } = s else { return Some(false) };
|
|
for q in [a, b] {
|
|
let point = plane.to_world(q) + plane.normal() * plane.height(from_p3(edge.front().point()));
|
|
let delta = point - origin;
|
|
let error =
|
|
1024. * f64::EPSILON * (1. + point.len() + origin.len()) + delta.x.abs() * ne.x + delta.y.abs() * ne.y + delta.z.abs() * ne.z;
|
|
if delta.dot(n).abs() + error > tol * 8. {
|
|
return Some(false);
|
|
}
|
|
}
|
|
}
|
|
return Some(true);
|
|
}
|
|
if let mt::Surface::RevolutedCurve(processor) = surface {
|
|
let r = processor.entity();
|
|
let mt::Curve::Line(line) = r.entity_curve() else { return None };
|
|
let o = from_p3(r.origin());
|
|
let axis = r.axis();
|
|
let axis = Vec3::new(axis.x, axis.y, axis.z).normalized()?;
|
|
let delta = from_p3(line.0) - o;
|
|
let radial = (delta - axis * delta.dot(axis)).len();
|
|
let u = axis.any_perp();
|
|
let v = axis.cross(u);
|
|
let matrix = *processor.transform();
|
|
let place = |p: Vec3, w: f64| {
|
|
let q = matrix * mt::Vector4::new(p.x, p.y, p.z, w);
|
|
(q.w == w).then_some(Vec3::new(q.x, q.y, q.z))
|
|
};
|
|
let center = plane.to_local(place(o, 1.)?);
|
|
let ha = plane.height(place(from_p3(line.0), 1.)?);
|
|
let hb = plane.height(place(from_p3(line.1), 1.)?);
|
|
if ha.min(hb) > levels.0 + tol * 8. || ha.max(hb) < levels.1 - tol * 8. {
|
|
return Some(false);
|
|
}
|
|
let (pu, pv) = (place(u, 0.)?, place(v, 0.)?);
|
|
let radius = pu.len() * radial;
|
|
if (pu.len() - pv.len()).abs() * radial > tol
|
|
|| pu.dot(pv).abs() * radial > tol
|
|
|| pu.dot(plane.normal()).abs() * radial > tol
|
|
|| pv.dot(plane.normal()).abs() * radial > tol
|
|
{
|
|
return None;
|
|
}
|
|
return Some(cap.iter().all(|s| matches!(*s,Seg2::Arc {center:c,radius:r,..} if c.dist(center)+(r-radius).abs()<=tol*8.)));
|
|
}
|
|
let (rows, cols) = match surface {
|
|
mt::Surface::BSplineSurface(s) => (s.control_points().len(), s.control_points().first()?.len()),
|
|
mt::Surface::NurbsSurface(s) => (s.control_points().len(), s.control_points().first()?.len()),
|
|
_ => return None,
|
|
};
|
|
*work = work.checked_sub(rows.checked_mul(cols)?.checked_mul(4)?)?;
|
|
let at = |i: usize, j: usize| -> Option<(Vec3, f64)> {
|
|
match surface {
|
|
mt::Surface::BSplineSurface(s) => s.control_points().get(i)?.get(j).map(|p| (from_p3(*p), 1.)),
|
|
mt::Surface::NurbsSurface(s) => s.control_points().get(i)?.get(j).map(|p| (Vec3::new(p.x / p.w, p.y / p.w, p.z / p.w), p.w)),
|
|
_ => None,
|
|
}
|
|
};
|
|
for transpose in [false, true] {
|
|
let point = |i, j| if transpose { at(j, i) } else { at(i, j) };
|
|
let (groups, length) = if transpose { (cols, rows) } else { (rows, cols) };
|
|
// The first measurement subset uses a degree-one axial ruling. Its
|
|
// positive rational interpolation is monotone and reaches both caps;
|
|
// projected incidence alone would miss an axially displaced support.
|
|
if !(2..=4096).contains(&groups) || length != 2 {
|
|
continue;
|
|
}
|
|
let axial_knots = match surface {
|
|
mt::Surface::BSplineSurface(s) => {
|
|
if transpose {
|
|
s.uknot_vec()
|
|
} else {
|
|
s.vknot_vec()
|
|
}
|
|
}
|
|
mt::Surface::NurbsSurface(s) => {
|
|
if transpose {
|
|
s.uknot_vec()
|
|
} else {
|
|
s.vknot_vec()
|
|
}
|
|
}
|
|
_ => return None,
|
|
};
|
|
if axial_knots.len() != 4
|
|
|| axial_knots.first() != axial_knots.get(1)
|
|
|| axial_knots.get(2) != axial_knots.last()
|
|
|| axial_knots.get(1)? >= axial_knots.get(2)?
|
|
{
|
|
continue;
|
|
}
|
|
let mut projected = Vec::new();
|
|
let mut valid = true;
|
|
let mut increasing = None;
|
|
for i in 0..groups {
|
|
let (base, w) = point(i, 0)?;
|
|
let h0 = plane.height(base);
|
|
let h1 = plane.height(point(i, 1)?.0);
|
|
if h0.min(h1) > levels.0 + tol * 8. || h0.max(h1) < levels.1 - tol * 8. || increasing.is_some_and(|s| s != (h1 > h0)) {
|
|
valid = false;
|
|
break;
|
|
}
|
|
increasing = Some(h1 > h0);
|
|
let p = plane.to_local(base);
|
|
for j in 0..length {
|
|
if plane.to_local(point(i, j)?.0).dist(p) > tol * 8. {
|
|
valid = false;
|
|
break;
|
|
}
|
|
}
|
|
let q = plane.to_world(p);
|
|
projected.push(mt::Vector4::new(q.x * w, q.y * w, q.z * w, w));
|
|
}
|
|
if !valid {
|
|
continue;
|
|
}
|
|
let knots = match surface {
|
|
mt::Surface::BSplineSurface(s) => {
|
|
if transpose {
|
|
s.vknot_vec()
|
|
} else {
|
|
s.uknot_vec()
|
|
}
|
|
}
|
|
mt::Surface::NurbsSurface(s) => {
|
|
if transpose {
|
|
s.vknot_vec()
|
|
} else {
|
|
s.uknot_vec()
|
|
}
|
|
}
|
|
_ => return None,
|
|
};
|
|
if knots.len() > 8192 {
|
|
return None;
|
|
}
|
|
*work = work.checked_sub(knots.len())?;
|
|
let c = mt::Curve::NurbsCurve(mt::NurbsCurve::new(mt::BSplineCurve::try_new(knots.clone(), projected).ok()?));
|
|
let (a, b) = c.range_tuple();
|
|
let support = mt::Edge::new(&mt::Vertex::new(c.subs(a)), &mt::Vertex::new(c.subs(b)), c);
|
|
let Some(fit) = crate::prism_curve::segment(&support, plane, tol, work) else {
|
|
continue;
|
|
};
|
|
let fit_loop = CapLoop { segments: fit.clone(), area: Interval::enclosing(0., 0.)? };
|
|
if let Some((center, radius)) = circle(&fit_loop, tol * 8.)
|
|
&& cap.iter().all(|s| matches!(*s,Seg2::Arc {center:c,radius:r,..} if c.dist(center)+(r-radius).abs()<=tol*8.))
|
|
{
|
|
return Some(true);
|
|
}
|
|
let same = loops_same(
|
|
&CapLoop { segments: cap.clone(), area: Interval::enclosing(0., 0.)? },
|
|
&CapLoop { segments: fit, area: Interval::enclosing(0., 0.)? },
|
|
tol * 8.,
|
|
work,
|
|
)?;
|
|
if same {
|
|
return Some(true);
|
|
}
|
|
}
|
|
Some(false)
|
|
}
|
|
|
|
fn circle(lp: &CapLoop, tolerance: f64) -> Option<(Vec2, f64)> {
|
|
let Seg2::Arc { center, radius, sweep: first, .. } = *lp.segments.first()? else { return None };
|
|
let mut sweep = 0.;
|
|
for s in &lp.segments {
|
|
let Seg2::Arc { center: c, radius: r, sweep: a, .. } = *s else { return None };
|
|
if c.dist(center) > tolerance || (r - radius).abs() > tolerance || a * first <= 0. {
|
|
return None;
|
|
}
|
|
sweep += a;
|
|
}
|
|
((sweep.abs() - std::f64::consts::TAU).abs() * radius <= tolerance).then_some((center, radius))
|
|
}
|
|
|
|
fn loops_same(a: &CapLoop, b: &CapLoop, tolerance: f64, work: &mut usize) -> Option<bool> {
|
|
*work = work.checked_sub(a.segments.len().checked_add(b.segments.len())?)?;
|
|
if (a.area.lo > 0.) != (b.area.lo > 0.) {
|
|
return Some(false);
|
|
}
|
|
if let (Some((ca, ra)), Some((cb, rb))) = (circle(a, tolerance), circle(b, tolerance)) {
|
|
return Some(ca.dist(cb) + (ra - rb).abs() <= tolerance);
|
|
}
|
|
if a.segments.len() != b.segments.len() {
|
|
return Some(false);
|
|
}
|
|
*work = work.checked_sub(a.segments.len().checked_mul(b.segments.len())?)?;
|
|
let mut used = vec![false; b.segments.len()];
|
|
for s in &a.segments {
|
|
let mut match_index = None;
|
|
for (i, t) in b.segments.iter().enumerate() {
|
|
if !used.get(i).copied()? && segment_same(s, t, tolerance) {
|
|
if match_index.is_some() {
|
|
return None;
|
|
}
|
|
match_index = Some(i);
|
|
}
|
|
}
|
|
let Some(i) = match_index else { return Some(false) };
|
|
*used.get_mut(i)? = true;
|
|
}
|
|
Some(true)
|
|
}
|
|
|
|
struct Certificate {
|
|
loops: Vec<CapLoop>,
|
|
volume: Interval,
|
|
lo: f64,
|
|
hi: f64,
|
|
}
|
|
|
|
fn preflight_body(body: &crate::Body) -> Option<()> {
|
|
if body.solid.boundaries().len() != 1 || body.mesh.is_some() {
|
|
return None;
|
|
}
|
|
let mut faces = 16_384usize;
|
|
let mut edges = 65_536usize;
|
|
for face in body.solid.face_iter() {
|
|
faces = faces.checked_sub(1)?;
|
|
// Borrow absolute wires: the oriented iterator would clone a hostile
|
|
// oversized wire before its edge count has been checked.
|
|
for wire in face.absolute_boundaries() {
|
|
for _ in wire.edge_iter() {
|
|
edges = edges.checked_sub(1)?;
|
|
}
|
|
}
|
|
}
|
|
Some(())
|
|
}
|
|
|
|
/// Initial material scope: a convex straight outer or a circle, with separated
|
|
/// circular holes. Wider profiles must acquire their own separation proof.
|
|
fn separated_material(loops: &[CapLoop], clearance: f64, work: &mut usize) -> Option<()> {
|
|
let outer = loops.iter().find(|l| l.area.lo > 0.)?;
|
|
if loops.iter().filter(|l| l.area.lo > 0.).count() != 1 {
|
|
return None;
|
|
}
|
|
let outer_circle = circle(outer, clearance);
|
|
if outer_circle.is_none() {
|
|
*work = work.checked_sub(outer.segments.len().checked_mul(outer.segments.len())?)?;
|
|
for s in &outer.segments {
|
|
let Seg2::Line { a, b } = *s else { return None };
|
|
let delta = b - a;
|
|
if delta.len() <= clearance {
|
|
return None;
|
|
}
|
|
for t in &outer.segments {
|
|
if delta.cross(t.start() - a) < -clearance * delta.len() {
|
|
return None;
|
|
}
|
|
}
|
|
}
|
|
let mut turn = 0.;
|
|
let mut angular_error = 0.;
|
|
for (i, s) in outer.segments.iter().enumerate() {
|
|
let Seg2::Line { a, b } = *s else {
|
|
return None;
|
|
};
|
|
let next = outer.segments.get((i + 1) % outer.segments.len())?;
|
|
let x = b - a;
|
|
let y = next.end() - next.start();
|
|
turn += x.cross(y).atan2(x.dot(y));
|
|
angular_error += clearance * 8. / x.len() + 256. * f64::EPSILON;
|
|
for (j, t) in outer.segments.iter().enumerate() {
|
|
if i != j && a.dist(t.start()) <= clearance {
|
|
return None;
|
|
}
|
|
}
|
|
}
|
|
if angular_error > 1e-3 || (turn - std::f64::consts::TAU).abs() > angular_error {
|
|
return None;
|
|
}
|
|
}
|
|
let mut holes = Vec::new();
|
|
for lp in loops {
|
|
if std::ptr::eq(lp, outer) {
|
|
continue;
|
|
}
|
|
if lp.area.hi >= 0. {
|
|
return None;
|
|
}
|
|
let (center, radius) = circle(lp, clearance)?;
|
|
if let Some((c, r)) = outer_circle {
|
|
if center.dist(c) + radius + clearance >= r {
|
|
return None;
|
|
}
|
|
} else {
|
|
*work = work.checked_sub(outer.segments.len())?;
|
|
for s in &outer.segments {
|
|
let Seg2::Line { a, b } = *s else { return None };
|
|
if (b - a).cross(center - a) / (b - a).len() <= radius + clearance {
|
|
return None;
|
|
}
|
|
}
|
|
}
|
|
*work = work.checked_sub(holes.len())?;
|
|
for (c, r) in &holes {
|
|
if center.dist(*c) <= radius + r + clearance {
|
|
return None;
|
|
}
|
|
}
|
|
holes.push((center, radius));
|
|
}
|
|
Some(())
|
|
}
|
|
|
|
fn certificate(body: &crate::Body, plane: &Plane, tol: f64) -> Option<Certificate> {
|
|
preflight_body(body)?;
|
|
let d = plane.normal();
|
|
let body_size = body.size();
|
|
let mut remaining = 16_384usize;
|
|
let mut caps = Vec::new();
|
|
let mut walls = Vec::new();
|
|
let mut edge_uses = 65_536usize;
|
|
for face in body.solid.face_iter() {
|
|
remaining = remaining.checked_sub(1)?;
|
|
for _ in face.edge_iter() {
|
|
edge_uses = edge_uses.checked_sub(1)?;
|
|
}
|
|
let surface = face.oriented_surface();
|
|
if let mt::Surface::Plane(p) = &surface {
|
|
let n = p.normal();
|
|
let n = Vec3::new(n.x, n.y, n.z).normalized()?;
|
|
if n.dot(d).abs() > 1. - 1e-12 {
|
|
if n.cross(d).len() * body_size > tol {
|
|
return None;
|
|
}
|
|
// A planar support's parameter origin may be arbitrarily far
|
|
// from its trimmed face. Use actual boundary geometry only.
|
|
let level = plane.height(from_p3(face.vertex_iter().next()?.point()));
|
|
caps.push((face, level, n.dot(d) > 0.));
|
|
if caps.len() > 2 {
|
|
return None;
|
|
}
|
|
continue;
|
|
}
|
|
}
|
|
walls.push((face, surface));
|
|
}
|
|
if caps.len() != 2 {
|
|
return None;
|
|
}
|
|
let top = caps.iter().find(|c| c.2)?;
|
|
let bottom = caps.iter().find(|c| !c.2)?;
|
|
let height = top.1 - bottom.1;
|
|
if height <= tol || !height.is_finite() {
|
|
return None;
|
|
}
|
|
for cap in [top, bottom] {
|
|
let mt::Surface::Plane(support) = cap.0.oriented_surface() else {
|
|
return None;
|
|
};
|
|
let origin = from_p3(support.origin());
|
|
let (n, ne) = plane_normal(&support)?;
|
|
for vertex in cap.0.vertex_iter() {
|
|
let point = from_p3(vertex.point());
|
|
let delta = point - origin;
|
|
let terms = (delta.x * n.x).abs() + (delta.y * n.y).abs() + (delta.z * n.z).abs();
|
|
let error = 1024. * f64::EPSILON * terms + delta.x.abs() * ne.x + delta.y.abs() * ne.y + delta.z.abs() * ne.z;
|
|
if delta.dot(n).abs() + error > tol * 2. || (plane.height(point) - cap.1).abs() > tol * 2. {
|
|
return None;
|
|
}
|
|
}
|
|
}
|
|
let mut work = 2_000_000;
|
|
let loops = top.0.boundaries().iter().map(|w| cap_loop(w, plane, tol, &mut work)).collect::<Option<Vec<_>>>()?;
|
|
let lower = bottom.0.boundaries().iter().map(|w| cap_loop(w, plane, tol, &mut work)).collect::<Option<Vec<_>>>()?;
|
|
if loops.len() != lower.len() || loops.is_empty() || loops.len() > 8192 {
|
|
return None;
|
|
}
|
|
let mut pair_work = 16_000_000;
|
|
for a in &loops {
|
|
// Bottom cap has the opposite material orientation.
|
|
let reverse = CapLoop { segments: a.segments.clone(), area: Interval { lo: -a.area.hi, hi: -a.area.lo } };
|
|
let mut matches = 0usize;
|
|
for b in &lower {
|
|
if loops_same(&reverse, b, tol * 8., &mut pair_work)? {
|
|
matches += 1;
|
|
}
|
|
}
|
|
if matches != 1 {
|
|
return None;
|
|
}
|
|
}
|
|
let mut controls = 262_144;
|
|
// Include cap fitting, support fitting and the <=8tol support-to-cap
|
|
// comparison, not just the control-net's projected ruling variation.
|
|
let mut corridor = tol * 16.;
|
|
let cap_edges = top.0.edge_iter().count().checked_add(bottom.0.edge_iter().count())?;
|
|
pair_work = pair_work.checked_sub((65_536 - edge_uses).checked_mul(cap_edges)?)?;
|
|
for (face, surface) in &walls {
|
|
let variation = wall_error(surface, d, height, &mut controls)?;
|
|
if variation > tol * 8. {
|
|
return None;
|
|
}
|
|
corridor = corridor.max(tol * 16. + variation);
|
|
if corridor > tol * 24. {
|
|
return None;
|
|
}
|
|
// Every wall must end on these actual cap wires. This is topology
|
|
// correspondence, not nearest geometry or stale face provenance.
|
|
let mut upper = 0usize;
|
|
let mut lower = 0usize;
|
|
for edge in face.edge_iter() {
|
|
if top.0.edge_iter().any(|e| e.id() == edge.id()) {
|
|
if !wall_supports_cap(surface, &edge, plane, (bottom.1, top.1), tol, &mut work)? {
|
|
return None;
|
|
}
|
|
upper += 1;
|
|
continue;
|
|
}
|
|
if bottom.0.edge_iter().any(|e| e.id() == edge.id()) {
|
|
if !wall_supports_cap(surface, &edge, plane, (bottom.1, top.1), tol, &mut work)? {
|
|
return None;
|
|
}
|
|
lower += 1;
|
|
continue;
|
|
}
|
|
let mt::Curve::Line(line) = edge.curve() else { return None };
|
|
if from_p3(line.0).dist(from_p3(edge.absolute_front().point())) > tol || from_p3(line.1).dist(from_p3(edge.absolute_back().point())) > tol
|
|
{
|
|
return None;
|
|
}
|
|
let delta = from_p3(line.1) - from_p3(line.0);
|
|
if delta.cross(d).len() > corridor || delta.dot(d).abs() < height - tol * 8. {
|
|
return None;
|
|
}
|
|
}
|
|
if upper == 0 || upper != lower {
|
|
return None;
|
|
}
|
|
}
|
|
let mut area = Interval::enclosing(0., 0.)?;
|
|
let mut length = 0.;
|
|
let mut count = 0usize;
|
|
for lp in &loops {
|
|
area = area.add(lp.area)?;
|
|
for edge in &lp.segments {
|
|
length += edge.length();
|
|
count = count.checked_add(1)?;
|
|
}
|
|
}
|
|
if count > 8192 || area.lo <= 0. {
|
|
return None;
|
|
}
|
|
separated_material(&loops, corridor * 16., &mut pair_work)?;
|
|
// Each boundary's monotone radial/linear corridor contributes at most a
|
|
// two-sided strip plus four endpoint disks. Unknown intersecting or
|
|
// touching corridors will be refused by the separation check below.
|
|
let strip = 2. * length * corridor + 4. * std::f64::consts::PI * count as f64 * corridor * corridor;
|
|
// Each cap's vertices are within 2tol and its whole rational edge
|
|
// control heights within another tol; include both cap corridors.
|
|
let h_error = tol * 8.;
|
|
if area.lo <= strip || height <= h_error {
|
|
return None;
|
|
}
|
|
let lo = (area.lo - strip) * (height - h_error);
|
|
let hi = (area.hi + strip) * (height + h_error);
|
|
let volume = Interval::enclosing((lo + hi) * 0.5, (hi - lo) * 0.5)?;
|
|
if volume.lo <= 0. || volume.hi - volume.lo > 0.001 {
|
|
return None;
|
|
}
|
|
Some(Certificate { loops, volume, lo: bottom.1, hi: top.1 })
|
|
}
|
|
|
|
/// Independently prove equality of the supported current uniform prisms.
|
|
/// `None` means unproved, including budget exhaustion; it never grants a copy.
|
|
pub fn certified_prism_equivalent(a: &crate::Body, b: &crate::Body) -> crate::Result<Option<bool>> {
|
|
crate::guard("certified prism comparison", || {
|
|
let result = (|| {
|
|
preflight_body(a)?;
|
|
preflight_body(b)?;
|
|
if !a.size().is_finite() || !b.size().is_finite() {
|
|
return None;
|
|
}
|
|
let size = a.size().max(b.size()).max(1.);
|
|
let tol = size * 1e-10;
|
|
let identity_tol = size * 1e-7;
|
|
let origin = a.solid.vertex_iter().next().map(|v| from_p3(v.point()))?;
|
|
let mut candidates = Vec::new();
|
|
let mut remaining = 16_384usize;
|
|
for f in a.solid.face_iter() {
|
|
remaining = remaining.checked_sub(1)?;
|
|
if let mt::Surface::Plane(p) = f.oriented_surface() {
|
|
let n = p.normal();
|
|
let n = Vec3::new(n.x, n.y, n.z).normalized()?;
|
|
if !candidates.iter().any(|d: &Vec3| d.dot(n).abs() > 1. - 1e-12) {
|
|
candidates.push(n);
|
|
if candidates.len() == 8 {
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
for d in candidates {
|
|
let plane = Plane::from_normal(origin, d)?;
|
|
let Some(ca) = certificate(a, &plane, tol) else { continue };
|
|
let Some(cb) = certificate(b, &plane, tol) else { continue };
|
|
if (ca.lo - cb.lo).abs() > identity_tol || (ca.hi - cb.hi).abs() > identity_tol || ca.loops.len() != cb.loops.len() {
|
|
return Some(false);
|
|
}
|
|
let mut work = 16_000_000;
|
|
let mut used = vec![false; cb.loops.len()];
|
|
for lp in &ca.loops {
|
|
let mut matched = None;
|
|
for (i, other) in cb.loops.iter().enumerate() {
|
|
if !used.get(i).copied()? && loops_same(lp, other, identity_tol, &mut work)? {
|
|
if matched.is_some() {
|
|
return None;
|
|
}
|
|
matched = Some(i);
|
|
}
|
|
}
|
|
let Some(i) = matched else { return Some(false) };
|
|
*used.get_mut(i)? = true;
|
|
}
|
|
let delta = (ca.volume.hi - cb.volume.lo).abs().max((cb.volume.hi - ca.volume.lo).abs());
|
|
return Some(delta <= 1e-5 * (1. + ca.volume.lo.abs()));
|
|
}
|
|
None
|
|
})();
|
|
Ok(result)
|
|
})
|
|
}
|
|
|
|
fn cap_loop(wire: &mt::Wire, plane: &Plane, tolerance: f64, work: &mut usize) -> Option<CapLoop> {
|
|
let mut segments = Vec::new();
|
|
let mut area = Interval::enclosing(0., 0.)?;
|
|
let mut ids = std::collections::HashSet::new();
|
|
for edge in wire.edge_iter() {
|
|
if !ids.insert(edge.id()) {
|
|
return None;
|
|
}
|
|
let certified = crate::prism_curve::segment(edge, plane, tolerance, work)?;
|
|
if segments.len().checked_add(certified.len())? > 8192 {
|
|
return None;
|
|
}
|
|
for segment in &certified {
|
|
area = area.add(green_interval(segment, tolerance)?)?;
|
|
}
|
|
segments.extend(certified);
|
|
}
|
|
if segments.is_empty() || !wire.is_closed() {
|
|
return None;
|
|
}
|
|
// Curve endpoints already agree with actual shared vertices in the
|
|
// whole-span proof. Check the fitted representation's continuity too.
|
|
for i in 0..segments.len() {
|
|
if segments.get(i)?.end().dist(segments.get((i + 1) % segments.len())?.start()) > 2. * tolerance {
|
|
return None;
|
|
}
|
|
}
|
|
Some(CapLoop { segments, area })
|
|
}
|
|
|
|
#[cfg(test)]
|
|
#[allow(clippy::unwrap_used, clippy::expect_used, clippy::panic)]
|
|
mod tests {
|
|
use super::*;
|
|
use mt::builder;
|
|
|
|
fn drilled(radius: f64, shift: f64) -> crate::Body {
|
|
use solidcraft_geom::{Loop2, Region2};
|
|
let region = Region2 {
|
|
outer: Loop2::polygon(&[Vec2::ZERO, Vec2::new(20., 0.), Vec2::new(20., 20.), Vec2::new(0., 20.)]),
|
|
holes: vec![Loop2::circle(Vec2::new(10. + shift, 10.), radius)],
|
|
};
|
|
crate::extrude(&Plane::XY, &[region], 0., 10.).unwrap().remove(0)
|
|
}
|
|
|
|
#[test]
|
|
fn actual_drilled_prism_and_step_reopen_are_certified_without_history() {
|
|
let source = drilled(6., 0.);
|
|
let before = crate::step_export(&[&source], "source").unwrap();
|
|
let imported = crate::step_import(&before).unwrap();
|
|
let copy = &imported.bodies[0].body;
|
|
let exact = 4000. - 360. * std::f64::consts::PI;
|
|
for body in [&source, copy] {
|
|
let cap = certificate(body, &Plane::XY, 3e-9).expect("actual complete current cap and walls");
|
|
assert!(cap.volume.lo <= exact && cap.volume.hi >= exact, "{:?}, expected {exact}", cap.volume);
|
|
assert!(cap.volume.hi - cap.volume.lo <= 0.001);
|
|
}
|
|
assert_eq!(certified_prism_equivalent(&source, copy).unwrap(), Some(true));
|
|
let wider = drilled(6.0001, 0.);
|
|
let a = certificate(&source, &Plane::XY, 3e-9).unwrap();
|
|
let b = certificate(&wider, &Plane::XY, 3e-9).unwrap();
|
|
assert!(a.volume.lo - b.volume.hi > 1e-5 * (1. + a.volume.hi));
|
|
assert_eq!(certified_prism_equivalent(&source, &wider).unwrap(), Some(false));
|
|
let shifted = drilled(6., 0.0001);
|
|
let same_volume = certificate(&shifted, &Plane::XY, 3e-9).unwrap();
|
|
assert!(same_volume.volume.lo <= exact && same_volume.volume.hi >= exact);
|
|
assert_eq!(certified_prism_equivalent(&source, &shifted).unwrap(), Some(false));
|
|
assert_eq!(crate::step_export(&[&source], "source").unwrap(), before);
|
|
}
|
|
|
|
#[test]
|
|
fn poorly_conditioned_plane_walls_do_not_get_certificates() {
|
|
let origin = mt::Point3::new(0., 0., 0.);
|
|
let u = mt::Point3::new(1., 0., 0.);
|
|
let v = mt::Point3::new(1., 0., 1e-9);
|
|
let surface = mt::Surface::Plane(mt::Plane::new(origin, u, v));
|
|
let edge = builder::line(&mt::Vertex::new(origin), &mt::Vertex::new(u));
|
|
let mut work = 2_000_000;
|
|
assert!(wall_error(&surface, Vec3::Z, 10., &mut work).is_none());
|
|
assert!(wall_supports_cap(&surface, &edge, &Plane::XY, (0., 10.), 1e-9, &mut work).is_none());
|
|
let normal = plane_normal(&mt::Plane::new(origin, u, mt::Point3::new(0., 0., 1.))).unwrap();
|
|
assert!(normal.0.dist(-Vec3::Y) < 1e-12);
|
|
assert!(normal.1.len() > 0.);
|
|
}
|
|
|
|
#[test]
|
|
fn rotated_current_boolean_prism_step_and_wrong_surface_incidence() {
|
|
let stock = crate::box_solid(Vec3::ZERO, Vec3::new(20., 20., 10.)).unwrap();
|
|
let cutter = crate::cylinder(Vec3::new(10., 10., 0.), Vec3::Z, 6., 10.).unwrap();
|
|
let source = crate::boolean(&stock, &cutter, crate::BoolOp::Cut).unwrap().unwrap();
|
|
let moved = crate::transform(&source, Vec3::new(30., 20., 10.), Vec3::ZERO, Vec3::new(1., 2., 3.), 0.7).unwrap();
|
|
let saved = crate::step_export(&[&moved], "actual boolean").unwrap();
|
|
let imported = crate::step_import(&saved).unwrap();
|
|
assert_eq!(certified_prism_equivalent(&moved, &imported.bodies[0].body).unwrap(), Some(true));
|
|
// Preserve every current cap/edge ID while replacing only the curved
|
|
// wall support with radius five. Topology incidence alone must fail.
|
|
let wrong = drilled(5., 0.);
|
|
let wrong_surface = wrong.solid.face_iter().find(|f| !matches!(f.surface(), mt::Surface::Plane(_))).unwrap().surface();
|
|
let faces = source
|
|
.solid
|
|
.face_iter()
|
|
.map(|f| {
|
|
if matches!(f.surface(), mt::Surface::Plane(_)) {
|
|
return f.clone();
|
|
}
|
|
let mut face = mt::Face::new(f.absolute_boundaries().clone(), wrong_surface.clone());
|
|
if !f.orientation() {
|
|
face.invert();
|
|
}
|
|
face
|
|
})
|
|
.collect::<Vec<_>>();
|
|
let invalid = crate::Body::new(crate::body::Solid::try_new(vec![faces.into()]).unwrap()).unwrap();
|
|
assert!(certificate(&invalid, &Plane::XY, 3e-9).is_none());
|
|
let actual = source.solid.face_iter().find(|f| !matches!(f.surface(), mt::Surface::Plane(_))).unwrap().surface();
|
|
let shifted = mt::Transformed::transformed(&actual, mt::Matrix4::from_translation(mt::Vector3::new(0., 0., 1000.)));
|
|
let top = source.solid.face_iter().find(|f| matches!(f.oriented_surface(),mt::Surface::Plane(p) if p.normal().z>0.9)).unwrap();
|
|
let arc = top.edge_iter().find(|e| !matches!(e.curve(), mt::Curve::Line(_))).unwrap();
|
|
assert_ne!(wall_supports_cap(&shifted, &arc, &Plane::XY, (0., 10.), 3e-9, &mut 2_000_000), Some(true));
|
|
}
|
|
|
|
#[test]
|
|
fn cap_support_origin_is_not_a_trimmed_face_height() {
|
|
let source = drilled(6., 0.);
|
|
let faces = source
|
|
.solid
|
|
.face_iter()
|
|
.map(|f| {
|
|
let mt::Surface::Plane(p) = f.surface() else {
|
|
return f.clone();
|
|
};
|
|
if p.normal().z.abs() < 0.9 {
|
|
return f.clone();
|
|
}
|
|
let origin = p.origin() + p.u_axis() * 1e8;
|
|
let plane = mt::Surface::Plane(mt::Plane::new(origin, origin + p.u_axis(), origin + p.v_axis()));
|
|
let mut face = mt::Face::new(f.absolute_boundaries().clone(), plane);
|
|
if !f.orientation() {
|
|
face.invert();
|
|
}
|
|
face
|
|
})
|
|
.collect::<Vec<_>>();
|
|
let same = crate::Body::new(crate::body::Solid::try_new(vec![faces.into()]).unwrap()).unwrap();
|
|
assert_eq!(certified_prism_equivalent(&source, &same).unwrap(), Some(true));
|
|
assert!(
|
|
loops_same(
|
|
&cap_loop(&source.solid.face_iter().last().unwrap().boundaries()[0], &Plane::XY, 3e-9, &mut 2_000_000).unwrap(),
|
|
&cap_loop(&source.solid.face_iter().last().unwrap().boundaries()[0], &Plane::XY, 3e-9, &mut 2_000_000).unwrap(),
|
|
1e-6,
|
|
&mut 0
|
|
)
|
|
.is_none()
|
|
);
|
|
}
|
|
|
|
#[test]
|
|
fn winding_retracing_and_uncertain_material_do_not_get_certificates() {
|
|
let arc = Seg2::Arc { center: Vec2::ZERO, radius: 6., start: 0., sweep: std::f64::consts::TAU };
|
|
let retraced = CapLoop { segments: vec![arc, arc.reversed(), arc], area: Interval::enclosing(36. * std::f64::consts::PI, 1e-8).unwrap() };
|
|
assert!(circle(&retraced, 1e-6).is_none());
|
|
let region = solidcraft_geom::Loop2::polygon(&[Vec2::ZERO, Vec2::new(20., 0.), Vec2::new(20., 20.), Vec2::new(0., 20.)]);
|
|
let doubled = CapLoop { segments: region.segs.iter().chain(®ion.segs).copied().collect(), area: Interval::enclosing(800., 1e-8).unwrap() };
|
|
assert!(separated_material(&[doubled], 1e-6, &mut 16_000_000).is_none());
|
|
let thin = crate::box_solid(Vec3::ZERO, Vec3::new(20., 20., 1e-5)).unwrap();
|
|
assert!(certificate(&thin, &Plane::XY, 3e-6).is_none());
|
|
}
|
|
|
|
#[test]
|
|
fn current_rational_circle_green_interval_encloses_independent_area() {
|
|
let r = 6.;
|
|
let vertices =
|
|
[mt::Point3::new(r, 0., 0.), mt::Point3::new(0., r, 0.), mt::Point3::new(-r, 0., 0.), mt::Point3::new(0., -r, 0.)].map(builder::vertex);
|
|
let mut edges = Vec::new();
|
|
for i in 0..4 {
|
|
let angle = (i as f64 + 0.5) * std::f64::consts::FRAC_PI_2;
|
|
edges.push(builder::circle_arc(&vertices[i], &vertices[(i + 1) % 4], mt::Point3::new(r * angle.cos(), r * angle.sin(), 0.)));
|
|
}
|
|
let wire: mt::Wire = edges.into();
|
|
let cap = cap_loop(&wire, &Plane::XY, 3e-9, &mut 2_000_000).unwrap();
|
|
let exact = std::f64::consts::PI * r * r;
|
|
assert!(cap.area.lo <= exact && cap.area.hi >= exact, "{:?}", cap.area);
|
|
assert!((cap.area.hi - cap.area.lo) * 10. < 0.001);
|
|
assert_eq!(cap.segments.len(), 4);
|
|
let reversed = cap_loop(&wire.inverse(), &Plane::XY, 3e-9, &mut 2_000_000).unwrap();
|
|
assert!(reversed.area.lo <= -exact && reversed.area.hi >= -exact);
|
|
}
|
|
}
|