Prepare 0.2.0 usability release and refresh desktop demo
This commit is contained in:
bherbruck committed
1 parent
a094fa46b5
commit
ab05b24dee
54 files changed
+4888
-234
No files matched your search
@@ -750,6 +750,7 @@ impl Body {
|
||||
return Ok(m);
|
||||
}
|
||||
m.seams = crate::topo::seam_flags(self, &m);
|
||||
(m.smooth_edges, m.sharp_edges) = crate::smooth_edges::flags(self, &m);
|
||||
Ok(m)
|
||||
}
|
||||
|
||||
|
||||
@@ -0,0 +1,186 @@
|
||||
//! Necessary boundary safety check for the narrowed-bore fallback.
|
||||
//! This is not a certificate of general Boolean equivalence. It prevents a
|
||||
//! cylinder recut from silently removing an additional material boundary that
|
||||
//! the fallback's coarse material samples do not see.
|
||||
use crate::body::{Body, from_p3};
|
||||
use solvecraft_geom::{Plane, Seg2, Vec3};
|
||||
use std::collections::HashSet;
|
||||
use truck_modeling as mt;
|
||||
|
||||
pub(super) fn boundaries_clear(body: &Body, bore_faces: &[usize], origin: Vec3, axis: Vec3, radius: f64, lo: f64, hi: f64, tol: f64) -> Option<()> {
|
||||
let plane = Plane::from_normal(origin, axis)?;
|
||||
if !origin.is_finite()
|
||||
|| !axis.is_finite()
|
||||
|| ![radius, lo, hi, tol].iter().all(|v| v.is_finite())
|
||||
|| radius <= tol
|
||||
|| hi <= lo
|
||||
|| tol <= 0.
|
||||
|| bore_faces.len() > 16_384
|
||||
{
|
||||
return None;
|
||||
}
|
||||
let mut remaining_faces = 16_384usize;
|
||||
let mut remaining_uses = 65_536usize;
|
||||
for face in body.solid.face_iter() {
|
||||
remaining_faces = remaining_faces.checked_sub(1)?;
|
||||
for wire in face.absolute_boundaries() {
|
||||
for _ in wire.edge_iter() {
|
||||
remaining_uses = remaining_uses.checked_sub(1)?;
|
||||
}
|
||||
}
|
||||
}
|
||||
// Exempt only existing coincident bore boundaries, never newly inferred
|
||||
// edges or every edge of a planar face adjoining a cylinder.
|
||||
let bore_faces: HashSet<usize> = bore_faces.iter().copied().collect();
|
||||
let mut bore_edges = HashSet::new();
|
||||
for (i, face) in body.solid.face_iter().enumerate() {
|
||||
if bore_faces.contains(&i) {
|
||||
for edge in face.edge_iter() {
|
||||
bore_edges.insert(edge.id());
|
||||
}
|
||||
}
|
||||
}
|
||||
let mut seen = HashSet::new();
|
||||
let mut work = 2_000_000usize;
|
||||
for edge in body.solid.edge_iter() {
|
||||
if !seen.insert(edge.id()) {
|
||||
continue;
|
||||
}
|
||||
let common_bore = bore_edges.contains(&edge.id());
|
||||
work = work.checked_sub(1)?;
|
||||
let curve = edge.curve();
|
||||
// Original matched-cylinder trims retain the legacy fallback policy.
|
||||
// Boolean-produced rational trims need not be exact analytic circles.
|
||||
// This exemption follows existing face/edge identity only; additional
|
||||
// trims must pass the whole-boundary proof below. It is deliberately
|
||||
// not a general proof that the cylinder contains no material.
|
||||
if common_bore && matches!(curve, mt::Curve::IntersectionCurve(_) | mt::Curve::BSplineCurve(_) | mt::Curve::NurbsCurve(_)) {
|
||||
continue;
|
||||
}
|
||||
if let mt::Curve::Line(line) = curve {
|
||||
let (a, b) = (from_p3(line.0), from_p3(line.1));
|
||||
if !a.is_finite() || !b.is_finite() {
|
||||
return None;
|
||||
}
|
||||
let ha = plane.height(a);
|
||||
let hb = plane.height(b);
|
||||
let error = tol + 512. * f64::EPSILON * (1. + a.len() + b.len() + origin.len());
|
||||
if common_bore {
|
||||
let pa = plane.to_local(a);
|
||||
let pb = plane.to_local(b);
|
||||
if pa.dist(pb) <= error && (pa.len() - radius).abs() + (pb.len() - radius).abs() <= error * 2. {
|
||||
continue;
|
||||
}
|
||||
return None;
|
||||
}
|
||||
let (mut t0, mut t1): (f64, f64) = (0., 1.);
|
||||
let dh = hb - ha;
|
||||
if dh.abs() <= error {
|
||||
if ha.min(hb) > hi + error || ha.max(hb) < lo - error {
|
||||
continue;
|
||||
}
|
||||
} else {
|
||||
let x = (lo - error - ha) / dh;
|
||||
let y = (hi + error - ha) / dh;
|
||||
t0 = t0.max(x.min(y));
|
||||
t1 = t1.min(x.max(y));
|
||||
if t0 > t1 {
|
||||
continue;
|
||||
}
|
||||
}
|
||||
let pa = plane.to_local(a + (b - a) * t0);
|
||||
let pb = plane.to_local(a + (b - a) * t1);
|
||||
let d = pb - pa;
|
||||
let t = if d.len2() > 0. { (-pa.dot(d) / d.len2()).clamp(0., 1.) } else { 0. };
|
||||
if (pa + d * t).len() <= radius + error {
|
||||
return None;
|
||||
}
|
||||
continue;
|
||||
}
|
||||
// Positive rational controls enclose the entire axial extent. Project
|
||||
// the current curve onto the axis plane; the existing whole-span
|
||||
// certificate proves its Line/Arc image, without point sampling.
|
||||
let (knots, controls): (mt::KnotVec, Vec<mt::Vector4>) = match curve {
|
||||
mt::Curve::BSplineCurve(c) => {
|
||||
if c.control_points().len() > 4096 || c.knot_vec().len() > 8192 {
|
||||
return None;
|
||||
}
|
||||
work = work.checked_sub(c.control_points().len().checked_add(c.knot_vec().len())?)?;
|
||||
(c.knot_vec().clone(), c.control_points().iter().map(|p| mt::Vector4::new(p.x, p.y, p.z, 1.)).collect())
|
||||
}
|
||||
mt::Curve::NurbsCurve(c) => {
|
||||
if c.control_points().len() > 4096 || c.knot_vec().len() > 8192 {
|
||||
return None;
|
||||
}
|
||||
work = work.checked_sub(c.control_points().len().checked_add(c.knot_vec().len())?)?;
|
||||
(c.knot_vec().clone(), c.control_points().clone())
|
||||
}
|
||||
_ => return None,
|
||||
};
|
||||
let mut heights = (f64::INFINITY, f64::NEG_INFINITY);
|
||||
let mut projection_error: f64 = 0.;
|
||||
let mut projected_origin = None;
|
||||
let mut projected_spread: f64 = 0.;
|
||||
let mut projected = Vec::with_capacity(controls.len());
|
||||
for p in controls {
|
||||
if ![p.x, p.y, p.z, p.w].iter().all(|v| v.is_finite()) || p.w <= 0. {
|
||||
return None;
|
||||
}
|
||||
let point = Vec3::new(p.x / p.w, p.y / p.w, p.z / p.w);
|
||||
if !point.is_finite() {
|
||||
return None;
|
||||
}
|
||||
let error = tol + 512. * f64::EPSILON * (1. + point.len() + origin.len());
|
||||
projection_error = projection_error.max(error);
|
||||
let h = plane.height(point);
|
||||
heights.0 = heights.0.min(h - error);
|
||||
heights.1 = heights.1.max(h + error);
|
||||
let q = plane.to_local(point);
|
||||
let first = *projected_origin.get_or_insert(q);
|
||||
projected_spread = projected_spread.max(q.dist(first));
|
||||
projected.push(mt::Vector4::new(q.x * p.w, q.y * p.w, 0., p.w));
|
||||
}
|
||||
if heights.0 > hi || heights.1 < lo {
|
||||
continue;
|
||||
}
|
||||
// A swept seam may be represented as a rational straight edge. Its
|
||||
// positive control hull proves a point-sized radial projection even
|
||||
// though a nondegenerate planar Line/Arc recognizer cannot accept it.
|
||||
if let Some(q) = projected_origin {
|
||||
if common_bore && projected_spread + projection_error + (q.len() - radius).abs() <= tol * 8. {
|
||||
continue;
|
||||
}
|
||||
if !common_bore && q.len() - projected_spread - projection_error > radius + tol * 4. {
|
||||
continue;
|
||||
}
|
||||
}
|
||||
use mt::{BoundedCurve, ParametricCurve};
|
||||
let projected = mt::Curve::NurbsCurve(mt::NurbsCurve::new(mt::BSplineCurve::try_new(knots, projected).ok()?));
|
||||
let (a, b) = projected.range_tuple();
|
||||
let projected_edge = mt::Edge::new(&mt::Vertex::new(projected.subs(a)), &mt::Vertex::new(projected.subs(b)), projected);
|
||||
let segments = crate::prism_curve::segment(&projected_edge, &Plane::XY, tol, &mut work)?;
|
||||
for segment in segments {
|
||||
if common_bore {
|
||||
if let Seg2::Arc { center, radius: r, .. } = segment
|
||||
&& center.len() + (r - radius).abs() + projection_error <= tol * 8.
|
||||
{
|
||||
continue;
|
||||
}
|
||||
return None;
|
||||
}
|
||||
let clearance = match segment {
|
||||
Seg2::Line { a, b } => {
|
||||
let d = b - a;
|
||||
let t = if d.len2() > 0. { (-a.dot(d) / d.len2()).clamp(0., 1.) } else { 0. };
|
||||
(a + d * t).len()
|
||||
}
|
||||
Seg2::Arc { center, radius, .. } => radius - center.len(),
|
||||
_ => return None,
|
||||
};
|
||||
if clearance <= radius + tol * 4. + projection_error {
|
||||
return None;
|
||||
}
|
||||
}
|
||||
}
|
||||
Some(())
|
||||
}
|
||||
@@ -0,0 +1,118 @@
|
||||
#![allow(clippy::unwrap_used, clippy::panic)]
|
||||
use crate::*;
|
||||
use solvecraft_geom::{Loop2, Plane, Region2, Seg2, Vec2, Vec3};
|
||||
|
||||
fn fixture() -> (Body, Body, Vec3, f64) {
|
||||
let slope = 1f64.to_radians().tan();
|
||||
let hs = |normal: Vec3, point: Vec3| {
|
||||
let n = normal.normalized().unwrap();
|
||||
HalfSpace { n, d: n.dot(point) }
|
||||
};
|
||||
let stock = convex_polyhedron(&[
|
||||
hs(Vec3::X, Vec3::new(40., 0., 0.)),
|
||||
hs(-Vec3::X, Vec3::ZERO),
|
||||
hs(Vec3::Y, Vec3::new(0., 30., 0.)),
|
||||
hs(-Vec3::Y, Vec3::ZERO),
|
||||
hs(-Vec3::Z, Vec3::ZERO),
|
||||
hs(Vec3::new(-slope, 0., 1.), Vec3::new(0., 0., 20.)),
|
||||
])
|
||||
.unwrap();
|
||||
let cutter = cylinder(Vec3::new(20., 15., -5.), Vec3::Z, 5., 40.).unwrap();
|
||||
let bored = boolean(&stock, &cutter, BoolOp::Cut).unwrap().unwrap();
|
||||
let web = box_solid(Vec3::new(14., 14., 16.05), Vec3::new(26., 16., 16.07)).unwrap();
|
||||
let target = boolean(&bored, &web, BoolOp::Union).unwrap().unwrap();
|
||||
let region = Region2 { outer: Loop2::circle(Vec2::new(20., 15.), 8.), holes: vec![Loop2::circle(Vec2::new(20., 15.), 5.).reversed()] };
|
||||
let ring = extrude(&Plane::XY, &[region], 14., 24.).unwrap().remove(0);
|
||||
let web_volume = 0.02 * 2. * (24f64.sqrt() + 25. * 0.2f64.asin());
|
||||
let fine_target = target.tessellate_with(1e-4, true).unwrap();
|
||||
let added = fine_target.measure().volume - bored.tessellate_with(1e-4, true).unwrap().measure().volume;
|
||||
assert!((added - web_volume).abs() < 0.005, "actual connected web adds {added}, expected {web_volume}");
|
||||
for x in [16., 20., 24.] {
|
||||
for y in [14.2, 15., 15.8] {
|
||||
assert!(fine_target.contains(Vec3::new(x, y, 16.06)), "source web must exist at {x},{y}");
|
||||
}
|
||||
}
|
||||
assert!(!fine_target.contains(Vec3::new(20., 12., 16.06)));
|
||||
let floor_volume = (20. + slope * 20.) * (1200. - 25. * std::f64::consts::PI);
|
||||
let ring_added = std::f64::consts::PI * (64. - 25.) * (24. - (20. + slope * 20.));
|
||||
(target, ring, Vec3::new(20., 15., 16.06), floor_volume + web_volume + ring_added)
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn thin_connected_web_between_bore_stations_survives_union() {
|
||||
let (target, ring, probe, expected) = fixture();
|
||||
let target_before = step_export(&[&target], "source").unwrap();
|
||||
let ring_before = step_export(&[&ring], "ring").unwrap();
|
||||
let size = target.size().max(ring.size());
|
||||
let tol = (size * 5e-4).max(size * 1e-5);
|
||||
let ma = target.tessellate_with(tol, true).unwrap();
|
||||
let mb = ring.tessellate_with(tol, true).unwrap();
|
||||
assert!(ma.contains(probe), "fixture must contain the real connected web");
|
||||
assert!(!mb.contains(probe));
|
||||
let va = ma.measure().volume;
|
||||
let vb = mb.measure().volume;
|
||||
let slack = 2e-3 * (va + vb) + 1e-9;
|
||||
let accept = |body: &Body| {
|
||||
let mr = body.tessellate_with(tol, true).unwrap();
|
||||
let v = mr.measure().volume;
|
||||
let mismatch = super::mismatch(&ma, &mb, &mr, BoolOp::Union);
|
||||
eprintln!("fallback candidate V={v}, expected={expected}, web={}, mismatch={mismatch:?}", mr.contains(probe));
|
||||
v > 0. && v >= va.max(vb) - slack && v <= va + vb + slack && mismatch.is_none_or(|f| f < 0.004)
|
||||
};
|
||||
let direct = crate::coplanar::bore_union(&target, &ring, &accept);
|
||||
eprintln!("direct fallback result={}", direct.is_some());
|
||||
let public = boolean(&target, &ring, BoolOp::Union);
|
||||
if let Ok(dir) = std::env::var("SOLVECRAFT_BORE_AUDIT_DIR") {
|
||||
let dir = std::path::Path::new(&dir);
|
||||
std::fs::create_dir_all(dir).unwrap();
|
||||
std::fs::write(dir.join("target-with-web.step"), step_export(&[&target], "target with web").unwrap()).unwrap();
|
||||
if let Some(result) = &direct {
|
||||
std::fs::write(dir.join("direct-union.step"), step_export(&[result], "direct union").unwrap()).unwrap();
|
||||
}
|
||||
if let Ok(Some(result)) = &public {
|
||||
std::fs::write(dir.join("public-union.step"), step_export(&[result], "public union").unwrap()).unwrap();
|
||||
}
|
||||
}
|
||||
let direct_keeps_web = direct.as_ref().is_none_or(|result| result.tessellate_with(1e-3, true).unwrap().contains(probe));
|
||||
let public_keeps_web = match &public {
|
||||
Ok(Some(result)) => {
|
||||
let keep = result.tessellate_with(1e-3, true).unwrap().contains(probe);
|
||||
eprintln!("public Union result web={keep}");
|
||||
keep
|
||||
}
|
||||
result => {
|
||||
eprintln!("public Union honestly refused: {result:?}");
|
||||
true
|
||||
}
|
||||
};
|
||||
assert!(direct_keeps_web, "accepted fallback removed real web material");
|
||||
assert!(public_keeps_web, "public Union removed real web material");
|
||||
assert_eq!(step_export(&[&target], "source").unwrap(), target_before);
|
||||
assert_eq!(step_export(&[&ring], "ring").unwrap(), ring_before);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn unproved_curved_boundary_between_clear_endpoints_refuses_recut() {
|
||||
let a = Vec2::new(6., -1.);
|
||||
let b = Vec2::new(6., 1.);
|
||||
let c = Vec2::new(8., 1.);
|
||||
let d = Vec2::new(8., -1.);
|
||||
let region = Region2 {
|
||||
outer: Loop2 {
|
||||
segs: vec![
|
||||
Seg2::Cubic { p0: a, p1: Vec2::new(3., -0.5), p2: Vec2::new(3., 0.5), p3: b },
|
||||
Seg2::Line { a: b, b: c },
|
||||
Seg2::Line { a: c, b: d },
|
||||
Seg2::Line { a: d, b: a },
|
||||
],
|
||||
}
|
||||
.reversed(),
|
||||
holes: vec![],
|
||||
};
|
||||
let body = extrude(&Plane::XY, &[region], 16.05, 16.07).unwrap().remove(0);
|
||||
assert!(a.len() > 5. && b.len() > 5., "endpoint samples miss the interior excursion");
|
||||
assert!(body.tessellate_with(1e-4, true).unwrap().contains(Vec3::new(4.5, 0., 16.06)));
|
||||
let before = step_export(&[&body], "curved source").unwrap();
|
||||
assert!(crate::bore_guard::boundaries_clear(&body, &[], Vec3::ZERO, Vec3::Z, 5., 14., 24., 1e-6).is_none());
|
||||
assert_eq!(step_export(&[&body], "curved source").unwrap(), before);
|
||||
}
|
||||
@@ -461,9 +461,11 @@ pub(crate) fn has_coincident_cylinders(a: &Body, b: &Body) -> bool {
|
||||
/// it is joined over). The mesh intersection can't take the coincident walls, so the tool's
|
||||
/// bore is narrowed (its walls pushed toward the axis, now clear of the target's), the bodies
|
||||
/// are joined, and the cylinder of the original bore is cut out again over the tool's length.
|
||||
/// That is exact when the cylinder holds no material of either body, which is checked on
|
||||
/// sample points. Each way (the tool's bore narrowed, else the target's) gives a result only
|
||||
/// when `accept` takes it. `None` when there are no coincident bores or no way works.
|
||||
/// Recutting is valid only where that cylinder holds no material of either body. The
|
||||
/// samples are plausibility checks; a whole-boundary guard rejects additional boundaries
|
||||
/// crossing the recut. Original common-bore trims retain the legacy eligibility policy.
|
||||
/// This is not a general emptiness or Boolean-equivalence certificate. Each direction
|
||||
/// returns a result only when `accept` takes it; otherwise this returns `None`.
|
||||
pub(crate) fn bore_union(a: &Body, b: &Body, accept: &dyn Fn(&Body) -> bool) -> Option<Body> {
|
||||
let size = a.size().max(b.size());
|
||||
let tol = (size * 1e-7).max(1e-9) * 10.0;
|
||||
@@ -536,6 +538,15 @@ pub(crate) fn bore_union(a: &Body, b: &Body, accept: &dyn Fn(&Body) -> bool) ->
|
||||
};
|
||||
let lo = if past(t0 - ext * 0.5) && past(t0 - ext) { t0 - ext } else { t0 };
|
||||
let hi = if past(t1 + ext * 0.5) && past(t1 + ext) { t1 + ext } else { t1 };
|
||||
// The samples above are plausibility checks, not proof that an
|
||||
// arbitrarily thin web is absent. Every additional current boundary
|
||||
// must be wholly clear of the recut; unsupported trims refuse.
|
||||
let other_bore: Vec<usize> = co.iter().filter(|x| same_cylinder(c, x, tol)).map(|x| x.index).collect();
|
||||
if crate::bore_guard::boundaries_clear(tool, &bore, c.origin, c.axis, c.radius, lo, hi, tol).is_none()
|
||||
|| crate::bore_guard::boundaries_clear(target, &other_bore, c.origin, c.axis, c.radius, lo, hi, tol).is_none()
|
||||
{
|
||||
continue;
|
||||
}
|
||||
let d = (c.radius * 0.2).min(size * 0.01);
|
||||
let run = || -> Result<Option<Body>> {
|
||||
let narrowed = crate::offset::offset_body(tool, |i, _| if bore.contains(&i) { d } else { 0.0 })?;
|
||||
|
||||
@@ -16,6 +16,7 @@ mod blend;
|
||||
mod blend_corner;
|
||||
mod blend_variable;
|
||||
mod body;
|
||||
mod bore_guard;
|
||||
mod build;
|
||||
mod coplanar;
|
||||
mod curveblend;
|
||||
@@ -39,8 +40,10 @@ mod polybool;
|
||||
mod polyhedron;
|
||||
mod prism;
|
||||
mod prism_curve;
|
||||
mod prism_measure;
|
||||
pub mod provenance;
|
||||
mod sew;
|
||||
mod smooth_edges;
|
||||
mod splitface;
|
||||
mod step;
|
||||
mod step_in;
|
||||
@@ -70,6 +73,7 @@ pub use ops::{BoolOp, boolean, can_join_touching_faces, split_by_plane, transfor
|
||||
pub use polybool::{faceted_boolean, planar_boolean};
|
||||
pub use polyhedron::{HalfSpace, convex_polyhedron, draft, offset_faces, shell};
|
||||
pub use prism::{extrude_planar_regions, planar_regions_boolean};
|
||||
pub use prism_measure::certified_prism_equivalent;
|
||||
pub use sew::{EdgeSpec, FaceSpec, SurfSpec, sew};
|
||||
pub use splitface::{SplitTool, split_body, split_faces, split_faces_with_map};
|
||||
pub use step::step_export;
|
||||
|
||||
@@ -7,6 +7,10 @@ use truck_modeling as mt;
|
||||
use crate::body::{Body, Solid, from_p3};
|
||||
use crate::{KernelError, Result, guard};
|
||||
|
||||
#[cfg(test)]
|
||||
#[path = "bore_union_tests.rs"]
|
||||
mod bore_union_tests;
|
||||
|
||||
#[derive(Clone, Copy, Debug, PartialEq, Eq, Serialize, Deserialize)]
|
||||
#[serde(rename_all = "snake_case")]
|
||||
pub enum BoolOp {
|
||||
|
||||
@@ -0,0 +1,903 @@
|
||||
//! Independent current-geometry intervals for uniform prisms.
|
||||
//! These are not mesh measurements or certificates of construction history.
|
||||
use crate::body::from_p3;
|
||||
use solvecraft_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 solvecraft_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 = solvecraft_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);
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,354 @@
|
||||
//! Display-only tangent joins. Unknown geometry retains its existing presentation.
|
||||
use solvecraft_geom::{Mesh, Vec3};
|
||||
use truck_modeling as mt;
|
||||
|
||||
use crate::body::{Body, from_p3};
|
||||
|
||||
const MAX_FACES: usize = 50_000;
|
||||
const MAX_EDGES: usize = 200_000;
|
||||
const WORK: usize = 2_000_000;
|
||||
const ANGLE: f64 = 1e-8;
|
||||
|
||||
#[derive(Clone, Copy)]
|
||||
struct Boundary {
|
||||
a: Vec3,
|
||||
b: Vec3,
|
||||
normal: Vec3,
|
||||
error: f64,
|
||||
}
|
||||
|
||||
fn charge(work: &mut usize, amount: usize) -> Option<()> {
|
||||
*work = work.checked_sub(amount)?;
|
||||
Some(())
|
||||
}
|
||||
|
||||
fn roundoff(points: &[Vec3]) -> f64 {
|
||||
points.iter().flat_map(|p| [p.x.abs(), p.y.abs(), p.z.abs()]).fold(1., f64::max) * f64::EPSILON * 64.
|
||||
}
|
||||
|
||||
/// Only a clamped Bezier profile swept by one common linear translation.
|
||||
/// Its whole ruling normal is bounded from the linear cross derivative, not sampled.
|
||||
fn ruled(surface: &mt::Surface, orientation: bool, work: &mut usize) -> Option<[Boundary; 2]> {
|
||||
let mt::Surface::BSplineSurface(surface) = surface else { return None };
|
||||
let controls = surface.control_points();
|
||||
let rows = controls.len();
|
||||
let cols = controls.first()?.len();
|
||||
if !(2..=33).contains(&rows) || !(2..=33).contains(&cols) || controls.iter().any(|row| row.len() != cols) {
|
||||
return None;
|
||||
}
|
||||
charge(work, rows.checked_mul(cols)?.checked_mul(4)?)?;
|
||||
let (ku, kv) = surface.knot_vecs();
|
||||
let clamped = |knots: &mt::KnotVec, count: usize| {
|
||||
let (Some(first), Some(last)) = (knots.first(), knots.last()) else { return false };
|
||||
knots.len() == count * 2
|
||||
&& first.is_finite()
|
||||
&& last.is_finite()
|
||||
&& first < last
|
||||
&& knots.iter().take(count).all(|k| k == first)
|
||||
&& knots.iter().skip(count).all(|k| k == last)
|
||||
};
|
||||
if !clamped(ku, rows) || !clamped(kv, cols) {
|
||||
return None;
|
||||
}
|
||||
for transpose in [false, true] {
|
||||
let (count, axial) = if transpose { (cols, rows) } else { (rows, cols) };
|
||||
if axial != 2 {
|
||||
continue;
|
||||
}
|
||||
let at = |i: usize, j: usize| {
|
||||
let p = if transpose { controls.get(j)?.get(i)? } else { controls.get(i)?.get(j)? };
|
||||
let p = from_p3(*p);
|
||||
p.is_finite().then_some(p)
|
||||
};
|
||||
let a = at(0, 0)?;
|
||||
let b = at(0, 1)?;
|
||||
let ruling = b - a;
|
||||
if ruling.len() < 1e-6 {
|
||||
continue;
|
||||
}
|
||||
let mut valid = true;
|
||||
for i in 0..count {
|
||||
let p = at(i, 0)?;
|
||||
let q = at(i, 1)?;
|
||||
if (q - p - ruling).len() > roundoff(&[a, b, p, q]) {
|
||||
valid = false;
|
||||
break;
|
||||
}
|
||||
}
|
||||
if !valid {
|
||||
continue;
|
||||
}
|
||||
let boundary = |index: usize, next: usize, end: bool| {
|
||||
let a = at(index, 0)?;
|
||||
let b = at(index, 1)?;
|
||||
let tangent = if end { a - at(next, 0)? } else { at(next, 0)? - a };
|
||||
let error = roundoff(&[a, b, at(next, 0)?]);
|
||||
// Ill-conditioned derivatives cannot certify a tiny angular join.
|
||||
if tangent.len() <= error / ANGLE || ruling.len() <= error / ANGLE {
|
||||
return None;
|
||||
}
|
||||
let cross = if transpose { ruling.cross(tangent) } else { tangent.cross(ruling) };
|
||||
let normal = cross.normalized()? * if orientation { 1. } else { -1. };
|
||||
let top_tangent = if end { b - at(next, 1)? } else { at(next, 1)? - b };
|
||||
let top_cross = if transpose { ruling.cross(top_tangent) } else { top_tangent.cross(ruling) };
|
||||
let top_normal = top_cross.normalized()? * if orientation { 1. } else { -1. };
|
||||
// Along this linear ruling the cross derivative is a linear
|
||||
// interpolation of these endpoint vectors. Positive alignment
|
||||
// excludes a singularity; all intermediate normal directions lie
|
||||
// between their directions. Include that whole-span variation.
|
||||
if normal.dot(top_normal) <= 0. {
|
||||
return None;
|
||||
}
|
||||
let conditioning = cross.len().min(top_cross.len());
|
||||
let normal_error = 8. * error * (tangent.len().max(top_tangent.len()) + ruling.len()) / conditioning + normal.dist(top_normal);
|
||||
if !normal_error.is_finite() || normal_error > ANGLE {
|
||||
return None;
|
||||
}
|
||||
Some(Boundary { a, b, normal, error: normal_error })
|
||||
};
|
||||
if let (Some(first), Some(last)) = (boundary(0, 1, false), boundary(count - 1, count - 2, true)) {
|
||||
return Some([first, last]);
|
||||
}
|
||||
}
|
||||
None
|
||||
}
|
||||
|
||||
fn on_boundary(boundaries: &[Boundary; 2], a: Vec3, b: Vec3) -> Option<Boundary> {
|
||||
boundaries.iter().find_map(|boundary| {
|
||||
let tolerance = roundoff(&[a, b, boundary.a, boundary.b]);
|
||||
((a.dist(boundary.a) <= tolerance && b.dist(boundary.b) <= tolerance) || (a.dist(boundary.b) <= tolerance && b.dist(boundary.a) <= tolerance))
|
||||
.then_some(*boundary)
|
||||
})
|
||||
}
|
||||
|
||||
pub(crate) fn flags(body: &Body, mesh: &Mesh) -> (Vec<bool>, Vec<bool>) {
|
||||
crate::guard("smooth display edges", || Ok(classify(body, mesh, WORK))).unwrap_or_default()
|
||||
}
|
||||
|
||||
fn classify(body: &Body, mesh: &Mesh, mut work: usize) -> (Vec<bool>, Vec<bool>) {
|
||||
if body.solid.face_iter().take(MAX_FACES + 1).count() > MAX_FACES || mesh.edges.len() > MAX_EDGES {
|
||||
return (Vec::new(), Vec::new());
|
||||
}
|
||||
let mut split_points = 0usize;
|
||||
for line in body.split_keep() {
|
||||
if charge(&mut work, line.len().saturating_add(1)).is_none() {
|
||||
return (Vec::new(), Vec::new());
|
||||
}
|
||||
let Some(total) = split_points.checked_add(line.len()) else { return (Vec::new(), Vec::new()) };
|
||||
split_points = total;
|
||||
}
|
||||
let mut boundaries = Vec::new();
|
||||
for face in body.solid.face_iter() {
|
||||
// Borrow the geometry through mapping, avoiding a deep clone of an
|
||||
// unknown/hostile surface before inspecting its control-net bounds.
|
||||
let wires = face.absolute_boundaries();
|
||||
if wires.len() > 16 {
|
||||
boundaries.push(None);
|
||||
continue;
|
||||
}
|
||||
let edge_count = wires.iter().fold(0usize, |sum, wire| sum.saturating_add(wire.len()));
|
||||
if edge_count > 16 || charge(&mut work, edge_count + 1).is_none() {
|
||||
boundaries.push(None);
|
||||
continue;
|
||||
}
|
||||
let mapped = face.try_mapped(|_| Some(()), |_| Some(()), |surface| ruled(surface, face.orientation(), &mut work));
|
||||
boundaries.push(mapped.map(|face| face.surface()));
|
||||
}
|
||||
let mut out = vec![false; mesh.edges.len()];
|
||||
let mut sharp = vec![false; mesh.edges.len()];
|
||||
// Same first-occurrence order as tessellate().edges; never remove/reindex it.
|
||||
let mut seen = std::collections::HashSet::new();
|
||||
let mut index = 0usize;
|
||||
for edge in body.solid.edge_iter() {
|
||||
if index >= mesh.edges.len() || seen.len() >= MAX_EDGES {
|
||||
break;
|
||||
}
|
||||
if charge(&mut work, 1).is_none() {
|
||||
break;
|
||||
}
|
||||
if !seen.insert(edge.id()) {
|
||||
continue;
|
||||
}
|
||||
let current = index;
|
||||
index += 1;
|
||||
let Some([a, b]) = mesh.edge_faces.get(current).map(Vec::as_slice) else { continue };
|
||||
if a == b {
|
||||
continue;
|
||||
}
|
||||
let Some(line) = mesh.edges.get(current) else { continue };
|
||||
let Some(split_work) = line.len().checked_add(1).and_then(|n| n.checked_mul(split_points)) else { break };
|
||||
if charge(&mut work, split_work).is_none() {
|
||||
break;
|
||||
}
|
||||
if crate::heal::on_split_line(line, body.split_keep()) {
|
||||
if let Some(flag) = sharp.get_mut(current) {
|
||||
*flag = true;
|
||||
}
|
||||
continue;
|
||||
}
|
||||
let (Some(Some(ba)), Some(Some(bb))) = (boundaries.get(*a as usize), boundaries.get(*b as usize)) else { continue };
|
||||
// Curve is inspected by borrowing. Only exact straight rulings qualify;
|
||||
// a bowed/twisted edge with matching endpoints must remain visible.
|
||||
let mapped = edge.try_mapped(
|
||||
|_| Some(()),
|
||||
|curve| match curve {
|
||||
mt::Curve::Line(line) => Some((from_p3(line.0), from_p3(line.1))),
|
||||
_ => None,
|
||||
},
|
||||
);
|
||||
let Some(mapped) = mapped else { continue };
|
||||
let (p, q) = mapped.curve();
|
||||
let (Some(na), Some(nb)) = (on_boundary(ba, p, q), on_boundary(bb, p, q)) else { continue };
|
||||
let uncertainty = na.error + nb.error;
|
||||
let cross = na.normal.cross(nb.normal).len();
|
||||
if na.normal.dot(nb.normal) > 0.
|
||||
&& cross + uncertainty <= ANGLE
|
||||
&& let Some(flag) = out.get_mut(current)
|
||||
{
|
||||
*flag = true;
|
||||
} else if (cross > ANGLE + uncertainty || na.normal.dot(nb.normal) < -uncertainty)
|
||||
&& let Some(flag) = sharp.get_mut(current)
|
||||
{
|
||||
*flag = true;
|
||||
}
|
||||
}
|
||||
(out, sharp)
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
#[allow(clippy::unwrap_used)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
use solvecraft_geom::{Loop2, Plane, Region2, Seg2, Vec2};
|
||||
|
||||
fn circle(radius: f64, count: usize) -> Loop2 {
|
||||
let step = std::f64::consts::TAU / count as f64;
|
||||
let k = 4. / 3. * (step / 4.).tan() * radius;
|
||||
let point = |angle: f64| Vec2::new(angle.cos(), angle.sin()) * radius;
|
||||
let tangent = |angle: f64| Vec2::new(-angle.sin(), angle.cos()) * k;
|
||||
Loop2 {
|
||||
segs: (0..count)
|
||||
.map(|i| {
|
||||
let a = i as f64 * step;
|
||||
let b = (i + 1) as f64 * step;
|
||||
let p0 = point(a);
|
||||
let p3 = if i + 1 == count { point(0.) } else { point(b) };
|
||||
Seg2::Cubic { p0, p1: p0 + tangent(a), p2: p3 - tangent(b), p3 }
|
||||
})
|
||||
.collect(),
|
||||
}
|
||||
}
|
||||
|
||||
fn ring(count: usize) -> Body {
|
||||
crate::extrude(&Plane::XY, &[Region2 { outer: circle(10., count), holes: vec![circle(9., count).reversed()] }], 0., 3.).unwrap().remove(0)
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn periodic_ring_hides_only_whole_smooth_wall_joins() {
|
||||
let body = ring(256);
|
||||
let mesh = body.display_mesh(0.02).unwrap();
|
||||
assert_eq!(mesh.edges.len(), 1536);
|
||||
assert_eq!(mesh.smooth_edges.iter().filter(|f| **f).count(), 512);
|
||||
for (index, (line, smooth)) in mesh.edges.iter().zip(&mesh.smooth_edges).enumerate() {
|
||||
let vertical = (line.first().unwrap().z - line.last().unwrap().z).abs() > 2.9;
|
||||
assert_eq!(*smooth, vertical);
|
||||
if !vertical {
|
||||
assert!(mesh.shaded_edge_visible(index));
|
||||
}
|
||||
}
|
||||
let raw = body.tessellate(0.02).unwrap();
|
||||
assert_eq!(mesh.edges, raw.edges);
|
||||
assert_eq!(mesh.edge_faces, raw.edge_faces);
|
||||
assert_eq!(mesh.tri_face, raw.tri_face);
|
||||
assert_eq!(mesh.triangles, raw.triangles);
|
||||
assert_eq!(mesh.positions, raw.positions);
|
||||
assert_eq!(mesh.normals, raw.normals);
|
||||
let moved = crate::transform_matrix(&body, [[1., 0., 0., 0.], [0., 0., 1., 0.], [0., -1., 0., 0.], [30., 20., 10., 1.]]).unwrap();
|
||||
assert_eq!(moved.display_mesh(0.02).unwrap().smooth_edges.iter().filter(|f| **f).count(), 512);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn authored_split_unknown_adjacency_and_exhausted_work_stay_visible() {
|
||||
let body = ring(8);
|
||||
let mesh = body.display_mesh(0.02).unwrap();
|
||||
let index = mesh.smooth_edges.iter().position(|f| *f).unwrap();
|
||||
let split = body.clone().with_split_lines(Some(std::sync::Arc::new(vec![mesh.edges[index].clone()])));
|
||||
assert!(!flags(&split, &mesh).0[index]);
|
||||
let mut unknown = mesh.clone();
|
||||
unknown.edge_faces[index].push(999);
|
||||
assert!(!flags(&body, &unknown).0[index]);
|
||||
assert!(classify(&body, &mesh, 0).0.iter().all(|f| !*f));
|
||||
// Matching endpoints alone cannot authorize a bowed current edge.
|
||||
let changed = Body::new(body.deep_copy()).unwrap();
|
||||
let edges = Body::unique_edges(&changed.solid);
|
||||
let edge = &edges[index];
|
||||
let (a, b) = (edge.front().point(), edge.back().point());
|
||||
let mid = mt::Point3::new((a.x + b.x) / 2. + 0.1, (a.y + b.y) / 2., (a.z + b.z) / 2.);
|
||||
edge.set_curve(mt::Curve::BSplineCurve(mt::BSplineCurve::new(mt::KnotVec::bezier_knot(2), vec![a, mid, b])));
|
||||
let (smooth, sharp) = flags(&changed, &mesh);
|
||||
assert!(!smooth[index] && !sharp[index]);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn a_small_real_crease_is_visible() {
|
||||
let mut outer = circle(10., 8);
|
||||
if let Seg2::Cubic { p0, p1, .. } = &mut outer.segs[0] {
|
||||
p1.x += (p1.y - p0.y) * 1e-6;
|
||||
}
|
||||
let body = crate::extrude(&Plane::XY, &[Region2 { outer, holes: vec![] }], 0., 3.).unwrap().remove(0);
|
||||
let mesh = body.display_mesh(0.02).unwrap();
|
||||
assert_eq!(mesh.smooth_edges.iter().filter(|f| **f).count(), 7);
|
||||
let crease = mesh
|
||||
.edges
|
||||
.iter()
|
||||
.position(|line| {
|
||||
line.first().is_some_and(|p| (p.x - 10.).abs() < 1e-8 && p.y.abs() < 1e-8)
|
||||
&& line.last().is_some_and(|p| (p.x - 10.).abs() < 1e-8 && p.y.abs() < 1e-8)
|
||||
&& line.first().zip(line.last()).is_some_and(|(a, b)| (a.z - b.z).abs() > 2.9)
|
||||
})
|
||||
.unwrap();
|
||||
assert!(!mesh.smooth_edges[crease]);
|
||||
assert!(mesh.sharp_edges[crease]);
|
||||
assert!(mesh.shaded_edge_visible(crease));
|
||||
assert!(mesh.edge_pickable(crease));
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn actual_authored_split_of_a_curved_ruled_wall_stays_visible() {
|
||||
let body = ring(4);
|
||||
let face = body
|
||||
.solid
|
||||
.face_iter()
|
||||
.position(|face| face.try_mapped(|_| Some(()), |_| Some(()), |surface| ruled(surface, face.orientation(), &mut { WORK })).is_some())
|
||||
.unwrap();
|
||||
let tool = crate::SplitTool::Plane(Plane::new(Vec3::new(7., 0., 0.), Vec3::Y, Vec3::Z).unwrap());
|
||||
let split = crate::split_faces(&body, &[face], &tool).unwrap();
|
||||
assert!(split.face_count() > body.face_count());
|
||||
assert!(!split.split_keep().is_empty());
|
||||
let mesh = split.display_mesh(0.02).unwrap();
|
||||
let authored: Vec<_> =
|
||||
mesh.edges.iter().enumerate().filter(|(_, line)| crate::heal::on_split_line(line, split.split_keep())).map(|(i, _)| i).collect();
|
||||
assert!(!authored.is_empty());
|
||||
for index in authored {
|
||||
assert!(!mesh.smooth_edges[index]);
|
||||
assert!(mesh.shaded_edge_visible(index));
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn singular_twisted_and_reversed_normal_certificates_are_conservative() {
|
||||
let surface = |points: Vec<Vec<mt::Point3>>| {
|
||||
mt::Surface::BSplineSurface(mt::BSplineSurface::new((mt::KnotVec::bezier_knot(points.len() - 1), mt::KnotVec::bezier_knot(1)), points))
|
||||
};
|
||||
let p = |x, y, z| mt::Point3::new(x, y, z);
|
||||
let straight = surface(vec![vec![p(0., 0., 0.), p(0., 0., 3.)], vec![p(1., 0., 0.), p(1., 0., 3.)]]);
|
||||
let forward = ruled(&straight, true, &mut { WORK }).unwrap();
|
||||
let backward = ruled(&straight, false, &mut { WORK }).unwrap();
|
||||
assert!(forward[0].normal.dot(backward[0].normal) < -0.99);
|
||||
let twisted = surface(vec![vec![p(0., 0., 0.), p(0., 0., 3.)], vec![p(1., 0., 0.), p(1., 1., 3.)]]);
|
||||
assert!(ruled(&twisted, true, &mut { WORK }).is_none());
|
||||
let singular = surface(vec![vec![p(0., 0., 0.), p(0., 0., 3.)], vec![p(0., 0., 0.), p(0., 0., 3.)], vec![p(1., 0., 0.), p(1., 0., 3.)]]);
|
||||
assert!(ruled(&singular, true, &mut { WORK }).is_none());
|
||||
assert!(ruled(&straight, true, &mut 0).is_none());
|
||||
}
|
||||
}
|
||||
Reference in new issue
Block a user