mirror of
https://github.com/bherbruck/solidcraft.git
synced 2026-11-03 10:35:58 +00:00
kernel: the mesher plans a schema 4 record, with no truck types
tess::plan takes a brep::Brep: faces, loops, shells' face sides and edges in the record's order. Surfaces go to the evaluator's form from the record (planes, splines and revolutions through SurfaceEval; a cylinder, cone, sphere or torus as the exact revolution of its line or rational-arc generator about its axis). Edges are sampled on the record's own curve evaluator (CurveEval), an intersection curve included, so tess/ssi.rs and the per-node parameters it carried are gone; every boundary node's parameters come from projection. No fallback is left in the planner's path: a curve that does not evaluate on its interval is an error naming the edge (it was a two-point chord, or truck's evaluation); a surface the evaluator cannot express an error naming the face and its kind. Errors name the face, its surface kind and, for a node off the surface, the loop, node and point. The module doc states the order: a cheaper layout where the face suits one, else constrained Delaunay in (u, v), else an error. MAX_ANGLE_BOUNDARY and MAX_RELAX are tess::size's own. Removed with truck: the winding test and pole-to-pole seam test that chose between this planner and truck's mesher, and boundaries_on_surfaces (truck blends). The tests in tess/ still build truck bodies and are left for the record-based Body. (not built: base in progress) tess/ has no errors under cargo check; the rest of the crate is the no-truck branch's. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_017itGh5fpMEousnZ91cYcVt
This commit is contained in:
jyviko
@@ -2,51 +2,37 @@
|
||||
//! same boundary points (watertight by construction).
|
||||
|
||||
use solidcraft_geom::Vec3;
|
||||
use truck_meshalgo::prelude::*;
|
||||
use truck_modeling as mt;
|
||||
|
||||
use solidcraft_geom::tess::EvalSurface;
|
||||
use solidcraft_geom::brep::eval::CurveEval;
|
||||
use solidcraft_geom::brep::{self, Curve};
|
||||
|
||||
use super::size::Rules;
|
||||
use crate::body::from_p3;
|
||||
use crate::{KernelError, Result};
|
||||
|
||||
/// Most nodes on one edge.
|
||||
pub(crate) const MAX_EDGE_NODES: usize = 20_000;
|
||||
/// Deepest bisection of one edge segment.
|
||||
const MAX_DEPTH: u32 = 24;
|
||||
|
||||
/// An edge's nodes, from its absolute front to its absolute back.
|
||||
/// An edge's nodes, from its start vertex (interval min) to its end vertex (interval max).
|
||||
pub(crate) struct Nodes {
|
||||
pub points: Vec<Vec3>,
|
||||
/// For an intersection edge: its two surfaces and, per node, both surfaces' parameters
|
||||
/// (u1, v1, u2, v2) where Newton placed the node (`None` where truck's evaluation did).
|
||||
pub on: Option<([EvalSurface; 2], Vec<Option<[f64; 4]>>)>,
|
||||
}
|
||||
|
||||
impl Nodes {
|
||||
fn plain(points: Vec<Vec3>) -> Nodes {
|
||||
Nodes { points, on: None }
|
||||
}
|
||||
|
||||
/// The parameters of each node on `surface`, where the edge knows them: `None` unless
|
||||
/// `surface` is one of the edge's two surfaces.
|
||||
pub(crate) fn params_on(&self, surface: &EvalSurface) -> Option<Vec<Option<[f64; 2]>>> {
|
||||
let (surfaces, params) = self.on.as_ref()?;
|
||||
let side = surfaces.iter().position(|s| s == surface)?;
|
||||
Some(params.iter().map(|q| q.map(|q| [q[2 * side], q[2 * side + 1]])).collect())
|
||||
}
|
||||
}
|
||||
|
||||
/// Nodes of `edge` from its absolute front to its absolute back: every chord within the
|
||||
/// tolerance of the curve (checked at the quarter points, so an S-bend in one segment is not
|
||||
/// missed), at most `h` long, and within the rules' facet angle. The end nodes are the
|
||||
/// vertices' own points, so edges meeting at a vertex meet exactly.
|
||||
pub(crate) fn sample_edge(edge: &mt::Edge, rules: &Rules, h: f64) -> Nodes {
|
||||
/// Nodes of edge `ei` of `b` from its start to its end: every chord within the tolerance of
|
||||
/// the curve (checked at the quarter points, so an S-bend in one segment is not missed), at
|
||||
/// most `h` long, and within the rules' facet angle. The end nodes are the vertices' own
|
||||
/// points, so edges meeting at a vertex meet exactly. A curve that does not evaluate somewhere
|
||||
/// on its interval is an error naming the edge.
|
||||
pub(crate) fn sample_edge(b: &brep::Brep, ei: usize, rules: &Rules, h: f64) -> Result<Nodes> {
|
||||
let fail = |why: String| KernelError::Failed(format!("tessellation: edge {ei}: {why}"));
|
||||
let edge = b.edges.get(ei).ok_or_else(|| fail("no such edge".into()))?;
|
||||
let vertex = |k: usize| {
|
||||
let v = edge.vertices.get(k).and_then(|&v| b.vertices.get(v as usize)).ok_or_else(|| fail("no such vertex".into()))?;
|
||||
Ok::<Vec3, KernelError>(Vec3::new(v.point[0], v.point[1], v.point[2]))
|
||||
};
|
||||
let (front, back) = (vertex(0)?, vertex(1)?);
|
||||
let tol = rules.tol;
|
||||
let front = from_p3(edge.absolute_front().point());
|
||||
let back = from_p3(edge.absolute_back().point());
|
||||
let curve = edge.curve();
|
||||
if matches!(curve, mt::Curve::Line(_)) {
|
||||
if matches!(edge.curve, Curve::Line { .. }) {
|
||||
// Evenly, as few segments as keep each at most `h`. The ends are the vertices' own
|
||||
// points: `front + (back − front) · 1` is not `back` to the last bit, and the face that
|
||||
// runs the edge backwards starts its loop with that node, where the face running it
|
||||
@@ -56,92 +42,65 @@ pub(crate) fn sample_edge(edge: &mt::Edge, rules: &Rules, h: f64) -> Nodes {
|
||||
if let Some(p) = pts.last_mut() {
|
||||
*p = back;
|
||||
}
|
||||
return Nodes::plain(pts);
|
||||
return Ok(Nodes { points: pts });
|
||||
}
|
||||
let (t0, t1) = curve.range_tuple();
|
||||
let [t0, t1] = edge.interval;
|
||||
if !(t0.is_finite() && t1.is_finite() && t1 > t0) {
|
||||
return Nodes::plain(vec![front, back]);
|
||||
return Err(fail(format!("interval [{t0}, {t1}]")));
|
||||
}
|
||||
// Our evaluators: the curve's own form where it has one (no allocation per point); for an
|
||||
// intersection curve, Newton on both surfaces from seeded parameters (see `ssi`), which
|
||||
// also gives each node's parameters on both faces; truck's evaluation only where that does
|
||||
// not converge.
|
||||
let fast = super::surface::eval_curve(&curve);
|
||||
let ssi = if fast.is_none() { super::ssi::Intersection::new(&curve) } else { None };
|
||||
type Pt = (Vec3, Option<[f64; 4]>);
|
||||
let at = |t: f64| -> Option<Pt> {
|
||||
match (&fast, &ssi) {
|
||||
(Some(c), _) => c.eval(t).map(|(p, _)| (p, None)),
|
||||
(None, Some(ix)) => ix.eval(t).map(|(p, q)| (p, Some(q))).or_else(|| Some((from_p3(curve.subs(t)), None))),
|
||||
(None, None) => Some((from_p3(curve.subs(t)), None)),
|
||||
}
|
||||
.filter(|(p, _)| p.is_finite())
|
||||
};
|
||||
let fallback = || Nodes::plain(vec![front, back]);
|
||||
let curve = CurveEval::new(&edge.curve).map_err(|e| fail(e.to_string()))?;
|
||||
let at =
|
||||
|t: f64| -> Result<Vec3> { curve.point(t).filter(|p| p.is_finite()).ok_or_else(|| fail(format!("the curve does not evaluate at t = {t}"))) };
|
||||
// A closed curve's chord between its ends is empty: start it from a few segments. An open
|
||||
// one starts whole, so an arc the tolerance leaves flat stays one segment.
|
||||
let initial: u32 = if edge.absolute_front().id() == edge.absolute_back().id() { 4 } else { 1 };
|
||||
let initial: u32 = if edge.vertices[0] == edge.vertices[1] { 4 } else { 1 };
|
||||
let t_at = |i: u32| t0 + (t1 - t0) * f64::from(i) / f64::from(initial);
|
||||
let mut ends: Vec<Pt> = Vec::with_capacity(initial as usize + 1);
|
||||
for i in 0..=initial {
|
||||
let Some(p) = at(t_at(i)) else { return fallback() };
|
||||
ends.push(p);
|
||||
}
|
||||
let ends: Vec<Vec3> = (0..=initial).map(|i| at(t_at(i))).collect::<Result<_>>()?;
|
||||
// Depth-first bisection, segment by segment, keeping the output ordered. Each segment
|
||||
// carries its end points and middle, so every parameter is evaluated once (an
|
||||
// intersection curve's point is a Newton solve).
|
||||
struct Seg {
|
||||
a: f64,
|
||||
b: f64,
|
||||
pa: Pt,
|
||||
pb: Pt,
|
||||
mid: Option<Pt>,
|
||||
pa: Vec3,
|
||||
pb: Vec3,
|
||||
mid: Option<Vec3>,
|
||||
depth: u32,
|
||||
}
|
||||
let Some(&first) = ends.first() else { return fallback() };
|
||||
let mut out: Vec<Pt> = vec![(front, first.1)];
|
||||
let mut out: Vec<Vec3> = vec![front];
|
||||
let mut stack: Vec<Seg> = Vec::with_capacity(64);
|
||||
for i in (0..initial).rev() {
|
||||
let (Some(pa), Some(pb)) = (<[Pt]>::get(&ends, i as usize).copied(), <[Pt]>::get(&ends, i as usize + 1).copied()) else {
|
||||
return fallback();
|
||||
};
|
||||
for (i, w) in ends.windows(2).enumerate().rev() {
|
||||
let (Some(&pa), Some(&pb)) = (w.first(), w.get(1)) else { continue };
|
||||
let i = i as u32;
|
||||
stack.push(Seg { a: t_at(i), b: t_at(i + 1), pa, pb, mid: None, depth: 0 });
|
||||
}
|
||||
while let Some(sg) = stack.pop() {
|
||||
let (a, b, pa, pb) = (sg.a, sg.b, sg.pa, sg.pb);
|
||||
let m = 0.5 * (a + b);
|
||||
let pm = match sg.mid {
|
||||
Some(p) => Some(p),
|
||||
None => at(m),
|
||||
Some(p) => p,
|
||||
None => at(m)?,
|
||||
};
|
||||
let (q1, q3) = (at(a + 0.25 * (b - a)), at(a + 0.75 * (b - a)));
|
||||
let (q1, q3) = (at(a + 0.25 * (b - a))?, at(a + 0.75 * (b - a))?);
|
||||
let done = sg.depth >= MAX_DEPTH || out.len() + stack.len() >= MAX_EDGE_NODES;
|
||||
let fine = || {
|
||||
let sag = [q1, pm, q3].iter().flatten().map(|q| q.0.dist_to_segment(pa.0, pb.0)).fold(0.0, f64::max);
|
||||
let short = !(pa.0.dist(pb.0) > h * (1.0 + 1e-9));
|
||||
let sag = [q1, pm, q3].iter().map(|q| q.dist_to_segment(pa, pb)).fold(0.0, f64::max);
|
||||
let short = !(pa.dist(pb) > h * (1.0 + 1e-9));
|
||||
// The two half chords of an arc turning θ differ by θ/2.
|
||||
let straight = match (rules.turn, pm) {
|
||||
(Some(t), Some(q)) => !t.exceeds(2.0 * turn(q.0 - pa.0, pb.0 - q.0), sag),
|
||||
_ => true,
|
||||
};
|
||||
let straight = rules.turn.is_none_or(|t| !t.exceeds(2.0 * turn(pm - pa, pb - pm), sag));
|
||||
short && straight && sag <= tol
|
||||
};
|
||||
match pm {
|
||||
Some(pm) if !done && !fine() => {
|
||||
stack.push(Seg { a: m, b, pa: pm, pb, mid: q3, depth: sg.depth + 1 });
|
||||
stack.push(Seg { a, b: m, pa, pb: pm, mid: q1, depth: sg.depth + 1 });
|
||||
}
|
||||
_ => out.push(pb),
|
||||
if !done && !fine() {
|
||||
stack.push(Seg { a: m, b, pa: pm, pb, mid: Some(q3), depth: sg.depth + 1 });
|
||||
stack.push(Seg { a, b: m, pa, pb: pm, mid: Some(q1), depth: sg.depth + 1 });
|
||||
} else {
|
||||
out.push(pb);
|
||||
}
|
||||
}
|
||||
if out.len() < 2 {
|
||||
return fallback();
|
||||
}
|
||||
if let Some(p) = out.last_mut() {
|
||||
p.0 = back;
|
||||
*p = back;
|
||||
}
|
||||
let on = ssi.map(|ix| (ix.surfaces().clone(), out.iter().map(|p| p.1).collect()));
|
||||
Nodes { points: out.into_iter().map(|p| p.0).collect(), on }
|
||||
Ok(Nodes { points: out })
|
||||
}
|
||||
|
||||
/// Angle (radians) between two directions; zero when either is empty.
|
||||
|
||||
+148
-382
@@ -1,20 +1,19 @@
|
||||
//! Display tessellation: each face's triangulation in its own parameters, with every edge
|
||||
//! sampled once so neighbouring faces share their boundary points exactly, handed over as a
|
||||
//! [`TessPlan`] for the CPU or GPU evaluator.
|
||||
//! SolidCraft's mesher: each face of a schema 4 record ([`brep::Brep`]) triangulated in its own
|
||||
//! parameters, with every edge sampled once so neighbouring faces share their boundary points
|
||||
//! exactly, handed over as a [`TessPlan`] for the CPU or GPU evaluator.
|
||||
//!
|
||||
//! Per face, cheapest first:
|
||||
//! - planar faces: a constrained Delaunay triangulation of the edge nodes alone (a plane needs
|
||||
//! no interior points);
|
||||
//! - faces whose boundary runs along their parameter box (cylinder halves, sphere and torus
|
||||
//! patches, fillets): a structured grid zipped to the edge nodes;
|
||||
//! - four-cornered faces whose sides are not parameter lines (mitred fillet strips, cross
|
||||
//! hole walls): a Coons blend of the four sides, zipped the same way;
|
||||
//! - anything else (trimmed faces, holes): constrained Delaunay with interior seeds and
|
||||
//! Ruppert refinement.
|
||||
//!
|
||||
//! A face the planner cannot handle (its surface has no evaluable form, a boundary that does
|
||||
//! not project or crosses itself) fails the whole plan, and the caller meshes the body the
|
||||
//! old way.
|
||||
//! Per face:
|
||||
//! 1. a cheaper layout where the face suits one: planar faces whose boundary runs along their
|
||||
//! parameter box, and curved ones, take a structured grid zipped to the edge nodes; a
|
||||
//! four-cornered face whose sides are not parameter lines (mitred fillet strips, cross hole
|
||||
//! walls) a Coons blend of its four sides; a trimmed face the patches of its cuts or the
|
||||
//! iso-line grid cut along its motorcycle traces;
|
||||
//! 2. otherwise, the general method for every face: constrained Delaunay in the face's (u, v),
|
||||
//! with interior seeds and Ruppert refinement, every node then evaluated on the exact
|
||||
//! surface (and refined where a triangle sags off it);
|
||||
//! 3. a face even that refuses (its loops cross or do not close in (u, v), or wind round a
|
||||
//! seam or pole) is an error naming the face and why. Its cause is upstream, in the face's
|
||||
//! boundary or parameters, and is fixed there; there is no other mesher to fall back to.
|
||||
|
||||
pub(crate) mod cdt;
|
||||
pub(crate) mod collar;
|
||||
@@ -28,17 +27,15 @@ pub(crate) mod pieces;
|
||||
pub(crate) mod sag;
|
||||
pub(crate) mod seam;
|
||||
pub(crate) mod size;
|
||||
pub(crate) mod ssi;
|
||||
pub(crate) mod surface;
|
||||
pub(crate) mod uv;
|
||||
|
||||
use std::collections::HashMap;
|
||||
|
||||
use solidcraft_geom::Vec3;
|
||||
use solidcraft_geom::brep::{self, Curve, Loop, Shell};
|
||||
use solidcraft_geom::tess::{EvalSurface, FacePlan, TessPlan};
|
||||
use truck_modeling as mt;
|
||||
|
||||
use crate::body::{Solid, edge_order};
|
||||
use crate::{KernelError, Result};
|
||||
|
||||
/// How the faces of a plan were meshed.
|
||||
@@ -81,63 +78,111 @@ fn unsupported(what: impl std::fmt::Display) -> KernelError {
|
||||
KernelError::Failed(format!("not supported yet: tessellation: {what}"))
|
||||
}
|
||||
|
||||
/// Plan the display tessellation of `solid` with chord tolerance `tol` and, when given, a
|
||||
/// facet angle no facet may turn more than (where it sags above the angle's floor).
|
||||
/// `boundary_tol`: how far the body's boundaries may lie off their surfaces (mm) as the file it
|
||||
/// was read from declares, `None` for a body this kernel built (see
|
||||
/// [`crate::body::Body::boundary_tol`]).
|
||||
pub(crate) fn plan(solid: &Solid, tol: f64, turn: Option<size::Turn>, boundary_tol: Option<f64>) -> Result<(TessPlan, PlanStats)> {
|
||||
let faces: Vec<&mt::Face> = solid.face_iter().collect();
|
||||
let (edges, edge_faces) = edge_order(faces.iter().copied());
|
||||
let rules = size::Rules::new(tol, size::extent(solid.vertex_iter().map(|v| crate::body::from_p3(v.point()))), turn, boundary_tol);
|
||||
/// The surface kind of a face, for an error that names it.
|
||||
fn kind(s: &brep::Surface) -> &'static str {
|
||||
match s {
|
||||
brep::Surface::Plane { .. } | brep::Surface::PlaneThrough { .. } => "plane",
|
||||
brep::Surface::Cylinder { .. } => "cylinder",
|
||||
brep::Surface::Cone { .. } => "cone",
|
||||
brep::Surface::Sphere { .. } => "sphere",
|
||||
brep::Surface::Torus { .. } => "torus",
|
||||
brep::Surface::Spline { .. } => "spline",
|
||||
brep::Surface::RationalSpline { .. } => "rational spline",
|
||||
brep::Surface::PlacedRevolution { .. } => "revolution",
|
||||
}
|
||||
}
|
||||
|
||||
/// Whether each face's normal is its surface's natural normal, as the shell that uses it says
|
||||
/// (a face no shell lists keeps its surface's).
|
||||
fn forward_faces(b: &brep::Brep) -> Vec<bool> {
|
||||
let mut forward = vec![true; b.faces.len()];
|
||||
for shell in &b.shells {
|
||||
if let Shell::Faces { faces, .. } = shell {
|
||||
for f in faces {
|
||||
if let Some(slot) = forward.get_mut(f.face as usize) {
|
||||
*slot = f.forward;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
forward
|
||||
}
|
||||
|
||||
/// The edges of each face's loops, in loop order.
|
||||
fn face_edges(f: &brep::Face) -> impl Iterator<Item = &brep::OrientedEdge> {
|
||||
f.loops.iter().flat_map(|l| match l {
|
||||
Loop::Edges { edges } => edges.as_slice(),
|
||||
Loop::Vertex { .. } => &[],
|
||||
})
|
||||
}
|
||||
|
||||
/// Plan the tessellation of `b` with chord tolerance `tol` and, when given, a facet angle no
|
||||
/// facet may turn more than (where it sags above the angle's floor). Faces and edges are the
|
||||
/// record's, in its order.
|
||||
pub(crate) fn plan(b: &brep::Brep, tol: f64, turn: Option<size::Turn>) -> Result<(TessPlan, PlanStats)> {
|
||||
let faces = &b.faces;
|
||||
let forward = forward_faces(b);
|
||||
// The faces each edge bounds, each once, in the order of the edge's uses.
|
||||
let edge_faces: Vec<Vec<usize>> = b
|
||||
.edges
|
||||
.iter()
|
||||
.map(|e| {
|
||||
let mut fs: Vec<usize> = Vec::with_capacity(e.uses.len());
|
||||
for u in &e.uses {
|
||||
if !fs.contains(&(u.face as usize)) {
|
||||
fs.push(u.face as usize);
|
||||
}
|
||||
}
|
||||
fs
|
||||
})
|
||||
.collect();
|
||||
let rules = size::Rules::new(tol, size::extent(b.vertices.iter().map(|v| Vec3::new(v.point[0], v.point[1], v.point[2]))), turn);
|
||||
// Each face's surface and target side length, then every edge once at the smaller size of
|
||||
// its faces; faces then read their nodes. A planar face costs microseconds and a curved
|
||||
// one about a hundred, so the pool is used only where curved faces (or intersection
|
||||
// curves) make up the work.
|
||||
let evals: Vec<Option<EvalSurface>> = faces.iter().map(|f| surface::eval_surface(&f.surface())).collect();
|
||||
let curved = evals.iter().filter(|e| !matches!(e, Some(EvalSurface::Plane { .. }))).count();
|
||||
let sizes =
|
||||
crate::par::map("face sizes", &evals, curved >= PARALLEL_SIZES, |eval| Ok(eval.as_ref().map(|s| size::face_size(s, &uv::domain(s), &rules))));
|
||||
let sized: Vec<(Option<EvalSurface>, Option<(f64, [f64; 2])>)> =
|
||||
evals.into_iter().zip(sizes).map(|(e, s)| s.map(|s| (e, s))).collect::<Result<_>>()?;
|
||||
let face_h = |f: usize| sized.get(f).and_then(|(_, s)| s.map(|(h, _)| h)).unwrap_or(rules.max_len);
|
||||
let evals: Vec<EvalSurface> = faces
|
||||
.iter()
|
||||
.enumerate()
|
||||
.map(|(fi, f)| surface::eval_surface(&f.surface, f.domain).map_err(|e| unsupported(format!("face {fi} ({}): {e}", kind(&f.surface)))))
|
||||
.collect::<Result<_>>()?;
|
||||
let curved = evals.iter().filter(|e| !matches!(e, EvalSurface::Plane { .. })).count();
|
||||
let sizes: Vec<(f64, [f64; 2])> =
|
||||
crate::par::map("face sizes", &evals, curved >= PARALLEL_SIZES, |s| Ok(size::face_size(s, &uv::domain(s), &rules)))
|
||||
.into_iter()
|
||||
.collect::<Result<_>>()?;
|
||||
let face_h = |f: usize| sizes.get(f).map_or(rules.max_len, |(h, _)| *h);
|
||||
let edge_h: Vec<f64> = edge_faces.iter().map(|fs| fs.iter().map(|&f| face_h(f)).fold(rules.max_len, f64::min)).collect();
|
||||
let jobs: Vec<(&mt::Edge, f64)> = edges.iter().zip(edge_h).collect();
|
||||
let solved = edges.iter().filter(|e| matches!(e.curve(), mt::Curve::IntersectionCurve(_))).count();
|
||||
let sampled = crate::par::map("edge nodes", &jobs, solved >= PARALLEL_CURVES || edges.len() >= PARALLEL_EDGES, |(e, h)| {
|
||||
Ok(edges::sample_edge(e, &rules, *h))
|
||||
});
|
||||
let mut nodes: HashMap<_, edges::Nodes> = HashMap::with_capacity(edges.len());
|
||||
for (e, n) in edges.iter().zip(sampled) {
|
||||
nodes.insert(e.id(), n?);
|
||||
}
|
||||
relax_crowded_faces(&faces, &jobs, &rules, &mut nodes);
|
||||
divide_collapsed_faces(&faces, &jobs, &rules, &mut nodes);
|
||||
let indexed: Vec<(usize, &mt::Face)> = faces.iter().copied().enumerate().collect();
|
||||
// A face's boundary nodes, loop by loop (each edge's last node is the next edge's first),
|
||||
// with their parameters on the face's surface where the edge knows them (an intersection
|
||||
// edge on this face's surface).
|
||||
type Boundary = (Vec<Vec<Vec3>>, Vec<Vec<Option<[f64; 2]>>>);
|
||||
let boundary = |face: &mt::Face, eval: &EvalSurface, nodes: &HashMap<mt::EdgeID, edges::Nodes>| -> Result<Boundary> {
|
||||
let jobs: Vec<(usize, f64)> = edge_h.into_iter().enumerate().collect();
|
||||
let solved = b.edges.iter().filter(|e| matches!(e.curve, Curve::SurfaceIntersection { .. })).count();
|
||||
let mut nodes: Vec<edges::Nodes> = crate::par::map("edge nodes", &jobs, solved >= PARALLEL_CURVES || jobs.len() >= PARALLEL_EDGES, |(e, h)| {
|
||||
edges::sample_edge(b, *e, &rules, *h)
|
||||
})
|
||||
.into_iter()
|
||||
.collect::<Result<_>>()?;
|
||||
relax_crowded_faces(b, &jobs, &rules, &mut nodes)?;
|
||||
let indexed: Vec<usize> = (0..faces.len()).collect();
|
||||
let planned = crate::par::map("face tessellation", &indexed, curved >= PARALLEL_FACES, |&fi| {
|
||||
let (Some(face), Some(eval), Some(size)) = (faces.get(fi), evals.get(fi), sizes.get(fi)) else {
|
||||
return Err(unsupported(format!("face {fi}: missing")));
|
||||
};
|
||||
// Boundary nodes, loop by loop (each edge's last node is the next edge's first). A
|
||||
// vertex loop (a pole) bounds no area and has no nodes of its own.
|
||||
let mut loops3: Vec<Vec<Vec3>> = Vec::new();
|
||||
let mut known: Vec<Vec<Option<[f64; 2]>>> = Vec::new();
|
||||
for wire in face.absolute_boundaries() {
|
||||
let (mut pts, mut uvs) = (Vec::new(), Vec::new());
|
||||
for e in wire.edge_iter() {
|
||||
let n = nodes.get(&e.id()).ok_or_else(|| unsupported("edge without nodes"))?;
|
||||
let (mut ep, mut eq) = (n.points.clone(), n.params_on(eval).unwrap_or_else(|| vec![None; n.points.len()]));
|
||||
if !e.orientation() {
|
||||
for l in &face.loops {
|
||||
let Loop::Edges { edges } = l else { continue };
|
||||
let mut pts = Vec::new();
|
||||
for oe in edges {
|
||||
let n = nodes.get(oe.edge as usize).ok_or_else(|| unsupported(format!("face {fi}: edge {} without nodes", oe.edge)))?;
|
||||
let mut ep = n.points.clone();
|
||||
if !oe.forward {
|
||||
ep.reverse();
|
||||
eq.reverse();
|
||||
}
|
||||
ep.pop();
|
||||
eq.pop();
|
||||
pts.extend(ep);
|
||||
uvs.extend(eq);
|
||||
}
|
||||
if pts.len() >= 2 {
|
||||
loops3.push(pts);
|
||||
known.push(uvs);
|
||||
}
|
||||
}
|
||||
Ok((loops3, known))
|
||||
@@ -157,95 +202,16 @@ pub(crate) fn plan(solid: &Solid, tol: f64, turn: Option<size::Turn>, boundary_t
|
||||
let (eval, size) = surface_of(fi)?;
|
||||
let (loops3, known) = boundary(face, &eval, nodes)?;
|
||||
let mut stats = PlanStats::default();
|
||||
let replan = |stats: &mut PlanStats| match plan_face(fi, &eval, &loops3, &known, rules, size, stats) {
|
||||
Ok(fp) => Ok(fp),
|
||||
// A boundary that doubles straight back on itself (a zero-width slit in the face,
|
||||
// which some STEP writers leave where two surfaces were trimmed apart) encloses no
|
||||
// area but crosses itself in (u, v). Dropping the slits leaves the same face.
|
||||
Err(first) => {
|
||||
let (slit3, slitq) = drop_slits(&loops3, &known, rules.tol * 1e-2);
|
||||
if slit3 == loops3 {
|
||||
return Err(first);
|
||||
}
|
||||
// Nothing left but slits: a sliver face (two edges between the same two vertices,
|
||||
// as STEP files write one) encloses no area, and the faces on either side already
|
||||
// meet along its nodes.
|
||||
if slit3.iter().all(|l| l.len() < 3) {
|
||||
return Ok(FacePlan { face: u32::try_from(fi).unwrap_or(u32::MAX), ..FacePlan::default() });
|
||||
}
|
||||
*stats = PlanStats::default();
|
||||
plan_face(fi, &eval, &slit3, &slitq, rules, size, stats).map_err(|_| first)
|
||||
}
|
||||
};
|
||||
let fp = replan(&mut stats)?;
|
||||
Ok((oriented_face(face, fp), eval, stats))
|
||||
};
|
||||
let mut planned = crate::par::map("face tessellation", &indexed, curved >= PARALLEL_FACES, |(fi, face)| plan_one(*fi, face, &nodes, &rules));
|
||||
// A face refused at this level of detail (a boundary that comes out degenerate, or that a
|
||||
// triangulation cannot follow, between nodes too far apart) has its edges divided again,
|
||||
// finer, and it and the faces beside those edges, which share their nodes, are planned
|
||||
// again from them.
|
||||
if planned.iter().any(Result::is_err) {
|
||||
let index_of: HashMap<mt::EdgeID, usize> = edges.iter().enumerate().map(|(i, e)| (e.id(), i)).collect();
|
||||
let face_edges: Vec<Vec<usize>> = faces
|
||||
.iter()
|
||||
.map(|f| f.absolute_boundaries().iter().flat_map(|w| w.edge_iter().filter_map(|e| index_of.get(&e.id()).copied())).collect())
|
||||
.collect();
|
||||
let mut finer: HashMap<usize, size::Rules> = HashMap::new();
|
||||
for factor in REFINE_STEPS {
|
||||
let refused: Vec<usize> = planned.iter().enumerate().filter(|(_, r)| r.is_err()).map(|(i, _)| i).collect();
|
||||
if refused.is_empty() {
|
||||
break;
|
||||
}
|
||||
let fine = rules.finer(factor);
|
||||
let mut touched: Vec<usize> = refused.iter().filter_map(|f| face_edges.get(*f)).flatten().copied().collect();
|
||||
touched.sort_unstable();
|
||||
touched.dedup();
|
||||
let resampled = crate::par::map("edge nodes", &touched, touched.len() >= PARALLEL_EDGES, |i| {
|
||||
Ok(jobs.get(*i).map(|(e, h)| (e.id(), edges::sample_edge(e, &fine, *h / factor))))
|
||||
});
|
||||
for (id, n) in resampled.into_iter().flatten().flatten() {
|
||||
nodes.insert(id, n);
|
||||
}
|
||||
finer.extend(refused.iter().map(|f| (*f, fine)));
|
||||
// The faces sharing a divided edge read its new nodes.
|
||||
let mut again: Vec<usize> = touched.iter().filter_map(|i| edge_faces.get(*i)).flatten().copied().collect();
|
||||
again.sort_unstable();
|
||||
again.dedup();
|
||||
let again: Vec<(usize, &mt::Face)> = again.into_iter().filter_map(|f| faces.get(f).map(|face| (f, *face))).collect();
|
||||
let results = crate::par::map("face tessellation", &again, again.len() >= PARALLEL_FACES, |(fi, face)| {
|
||||
plan_one(*fi, face, &nodes, finer.get(fi).unwrap_or(&rules))
|
||||
});
|
||||
for ((f, _), r) in again.iter().zip(results) {
|
||||
if let Some(slot) = planned.get_mut(*f) {
|
||||
*slot = r;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
// A face refused however finely its edges are divided is filled flat from its boundary
|
||||
// nodes (see `fallback`), so the rest of the body keeps its plan.
|
||||
for (fi, slot) in planned.iter_mut().enumerate() {
|
||||
if slot.is_ok() {
|
||||
continue;
|
||||
}
|
||||
let (Some(face), Ok((eval, _))) = (faces.get(fi), surface_of(fi)) else { continue };
|
||||
let Ok((loops3, _)) = boundary(face, &eval, &nodes) else { continue };
|
||||
if let Some(fp) = fallback::fill(fi, &eval, &loops3, off_surface(&rules), rules.tol) {
|
||||
*slot = Ok((oriented_face(face, fp), eval, PlanStats { filled: 1, ..PlanStats::default() }));
|
||||
}
|
||||
}
|
||||
let oriented = |e: &mt::Edge| -> Option<Vec<Vec3>> {
|
||||
let mut pts = nodes.get(&e.id())?.points.clone();
|
||||
if !e.orientation() {
|
||||
pts.reverse();
|
||||
}
|
||||
Some(pts)
|
||||
};
|
||||
let fp = plan_face(fi, kind(&face.surface), eval, &loops3, &rules, *size, &mut stats)?;
|
||||
let fwd = forward.get(fi).copied().unwrap_or(true);
|
||||
let fp =
|
||||
FacePlan { flip: !fwd, triangles: if fwd { fp.triangles } else { fp.triangles.into_iter().map(|[a, b, c]| [a, c, b]).collect() }, ..fp };
|
||||
Ok((fp, stats))
|
||||
});
|
||||
let mut plan = TessPlan::default();
|
||||
let mut stats = PlanStats::default();
|
||||
for r in planned {
|
||||
let (fp, eval, s) = r?;
|
||||
for (r, eval) in planned.into_iter().zip(evals) {
|
||||
let (fp, s) = r?;
|
||||
let surface = u32::try_from(plan.surfaces.len()).map_err(|_| unsupported("too many faces"))?;
|
||||
plan.surfaces.push(eval);
|
||||
plan.faces.push(FacePlan { surface: Some(surface), ..fp });
|
||||
@@ -254,239 +220,48 @@ pub(crate) fn plan(solid: &Solid, tol: f64, turn: Option<size::Turn>, boundary_t
|
||||
stats.delaunay += s.delaunay;
|
||||
stats.filled += s.filled;
|
||||
}
|
||||
plan.edges = edges.iter().filter_map(oriented).collect();
|
||||
plan.edges = nodes.into_iter().map(|n| n.points).collect();
|
||||
plan.edge_faces = edge_faces.into_iter().map(|fs| fs.into_iter().map(|f| u32::try_from(f).unwrap_or(u32::MAX)).collect()).collect();
|
||||
Ok((plan, stats))
|
||||
}
|
||||
|
||||
/// Whether every curved face's boundary lies on the face's surface, within `accept`: a point
|
||||
/// off the surface's own parameter range is off it. The planner needs this of every face; a
|
||||
/// blend that runs on past the end of a neighbouring round's surface does not have it.
|
||||
pub(crate) fn boundaries_on_surfaces(solid: &Solid, accept: f64) -> bool {
|
||||
use mt::{BoundedCurve, ParametricCurve};
|
||||
for face in solid.face_iter() {
|
||||
let Some(eval) = surface::eval_surface(&face.surface()) else { continue };
|
||||
if matches!(eval, EvalSurface::Plane { .. }) {
|
||||
continue;
|
||||
}
|
||||
let dom = uv::domain(&eval);
|
||||
let mut prev: Option<[f64; 2]> = None;
|
||||
for e in face.edge_iter() {
|
||||
let c = e.oriented_curve();
|
||||
let (lo, hi) = c.range_tuple();
|
||||
for k in 0..=32 {
|
||||
let p = crate::body::from_p3(c.subs(lo + (hi - lo) * f64::from(k) / 32.0));
|
||||
match uv::project(&eval, &dom, p, prev, accept) {
|
||||
Some(uv) => prev = Some(uv),
|
||||
None => return false,
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
true
|
||||
}
|
||||
|
||||
/// Bounded work for the facet angle, as the reference mesher bounds it: a face whose boundary
|
||||
/// would get more than `MAX_ANGLE_BOUNDARY` segments has its curved edges sampled again at twice
|
||||
/// the angle, as often as it takes, down to the tolerance alone (a plate with a thousand small
|
||||
/// holes is not worth a hundred times the work to draw each hole round). An edge keeps one set
|
||||
/// of nodes for both its faces, at the coarser of their levels, so every rim of a crowded face
|
||||
/// gets the same share whatever the order of the faces.
|
||||
fn relax_crowded_faces(faces: &[&mt::Face], jobs: &[(&mt::Edge, f64)], rules: &size::Rules, nodes: &mut HashMap<mt::EdgeID, edges::Nodes>) {
|
||||
use truck_meshalgo::tessellation::Deflection;
|
||||
let Some(turn) = rules.turn else { return };
|
||||
let h_of: HashMap<mt::EdgeID, (&mt::Edge, f64)> = jobs.iter().map(|(e, h)| (e.id(), (*e, *h))).collect();
|
||||
let boundaries: Vec<Vec<mt::EdgeID>> =
|
||||
faces.iter().map(|f| f.absolute_boundaries().iter().flat_map(|w| w.edge_iter().map(|e| e.id())).collect()).collect();
|
||||
let mut level: HashMap<mt::EdgeID, u8> = HashMap::new();
|
||||
for l in 1..=Deflection::MAX_RELAX {
|
||||
let segments = |id: &mt::EdgeID| nodes.get(id).map_or(0, |n| n.points.len().saturating_sub(1));
|
||||
let mut crowded: Vec<mt::EdgeID> = boundaries
|
||||
/// Bounded work for the facet angle: a face whose boundary would get more than
|
||||
/// [`size::MAX_ANGLE_BOUNDARY`] segments has its curved edges sampled again at twice the angle,
|
||||
/// as often as it takes, down to the tolerance alone (a plate with a thousand small holes is
|
||||
/// not worth a hundred times the work to draw each hole round). An edge keeps one set of nodes
|
||||
/// for both its faces, at the coarser of their levels, so every rim of a crowded face gets the
|
||||
/// same share whatever the order of the faces.
|
||||
fn relax_crowded_faces(b: &brep::Brep, jobs: &[(usize, f64)], rules: &size::Rules, nodes: &mut [edges::Nodes]) -> Result<()> {
|
||||
let Some(turn) = rules.turn else { return Ok(()) };
|
||||
let boundaries: Vec<Vec<usize>> = b.faces.iter().map(|f| face_edges(f).map(|oe| oe.edge as usize).collect()).collect();
|
||||
let curved = |e: usize| b.edges.get(e).is_some_and(|e| !matches!(e.curve, Curve::Line { .. }));
|
||||
let mut level: Vec<u8> = vec![0; nodes.len()];
|
||||
for l in 1..=size::MAX_RELAX {
|
||||
let segments = |e: &usize| nodes.get(*e).map_or(0, |n| n.points.len().saturating_sub(1));
|
||||
let mut crowded: Vec<usize> = boundaries
|
||||
.iter()
|
||||
.filter(|ids| ids.iter().map(segments).sum::<usize>() > Deflection::MAX_ANGLE_BOUNDARY)
|
||||
.filter(|es| es.iter().map(segments).sum::<usize>() > size::MAX_ANGLE_BOUNDARY)
|
||||
.flatten()
|
||||
.filter(|id| level.get(*id).copied().unwrap_or(0) < l)
|
||||
.filter(|id| h_of.get(*id).is_some_and(|(e, _)| !matches!(e.curve(), mt::Curve::Line(_))))
|
||||
.filter(|e| level.get(**e).copied().unwrap_or(0) < l && curved(**e))
|
||||
.copied()
|
||||
.collect();
|
||||
crowded.sort_unstable_by_key(|id| format!("{id:?}"));
|
||||
crowded.sort_unstable();
|
||||
crowded.dedup();
|
||||
if crowded.is_empty() {
|
||||
break;
|
||||
}
|
||||
let relaxed = size::Rules { turn: turn.relaxed(l), ..*rules };
|
||||
for id in crowded {
|
||||
if let Some((e, h)) = h_of.get(&id) {
|
||||
nodes.insert(id, edges::sample_edge(e, &relaxed, *h));
|
||||
level.insert(id, l);
|
||||
for e in crowded {
|
||||
let Some(&(_, h)) = jobs.get(e) else { continue };
|
||||
if let Some(slot) = nodes.get_mut(e) {
|
||||
*slot = edges::sample_edge(b, e, &relaxed, h)?;
|
||||
}
|
||||
if let Some(slot) = level.get_mut(e) {
|
||||
*slot = l;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Every loop with its zero-width slits removed: wherever the boundary steps out to a node and
|
||||
/// straight back to where it came from, the pair is dropped. A slit encloses no area, so the
|
||||
/// face is the same without it, and without it the boundary no longer runs along itself. The
|
||||
/// two sides need only agree within `eps` (a hundredth of the chord tolerance): a STEP writer
|
||||
/// may draw them as two curves a few nanometres apart, which no mesh could tell apart anyway.
|
||||
fn drop_slits(loops3: &[Vec<Vec3>], known: &[Vec<Option<[f64; 2]>>], eps: f64) -> (Vec<Vec<Vec3>>, Vec<Vec<Option<[f64; 2]>>>) {
|
||||
/// Slits dropped from one loop before giving up (the search is a scan per slit).
|
||||
const MAX_SLITS: usize = 4096;
|
||||
let mut out3: Vec<Vec<Vec3>> = Vec::with_capacity(loops3.len());
|
||||
let mut outq: Vec<Vec<Option<[f64; 2]>>> = Vec::with_capacity(loops3.len());
|
||||
for (li, pts) in loops3.iter().enumerate() {
|
||||
let mut pts = pts.clone();
|
||||
let mut qs = known.get(li).cloned().unwrap_or_else(|| vec![None; pts.len()]);
|
||||
qs.resize(pts.len(), None);
|
||||
for _ in 0..MAX_SLITS {
|
||||
let n = pts.len();
|
||||
if n < 3 {
|
||||
break;
|
||||
}
|
||||
let same = |a: usize, b: usize| pts.get(a).zip(pts.get(b)).is_some_and(|(p, q)| p.dist(*q) <= eps);
|
||||
let tip = (0..n).find(|&k| same((k + n - 1) % n, (k + 1) % n));
|
||||
let Some(k) = tip else { break };
|
||||
// Drop the tip and the node that repeats the one before it.
|
||||
let mut drop = [k, (k + 1) % n];
|
||||
drop.sort_unstable();
|
||||
for &i in drop.iter().rev() {
|
||||
if i < pts.len() {
|
||||
pts.remove(i);
|
||||
}
|
||||
if i < qs.len() {
|
||||
qs.remove(i);
|
||||
}
|
||||
}
|
||||
}
|
||||
if pts.len() >= 2 {
|
||||
out3.push(pts);
|
||||
outq.push(qs);
|
||||
}
|
||||
}
|
||||
(out3, outq)
|
||||
}
|
||||
|
||||
/// Boundary nodes a loop keeps only one of: nodes whose parameters land on the node before them.
|
||||
///
|
||||
/// A boundary that runs a little off its surface (a STEP file writes its curves and surfaces to
|
||||
/// its own uncertainty, so an edge may lie micrometres off the face it bounds) projects onto the
|
||||
/// nearest surface point. Around a convex corner of the surface every point of such an edge is
|
||||
/// nearest that corner, so a whole run of nodes comes out with the corner's parameters. The finer
|
||||
/// the chord tolerance, the more nodes the run holds, and a boundary that repeats a point in
|
||||
/// (u, v) encloses no area there and cannot be triangulated (#1721).
|
||||
///
|
||||
/// One node stands for the run. The nodes dropped are nearer to it than `SAME_NODE` on the
|
||||
/// surface, so the face is the same without them.
|
||||
fn drop_collapsed(l: &mut seam::Loop, speed: [f64; 2]) {
|
||||
let eps = |d: usize| {
|
||||
let s = speed.get(d).copied().unwrap_or(1.0);
|
||||
if s.is_finite() && s > 0.0 { SAME_NODE / s } else { 0.0 }
|
||||
};
|
||||
let (eu, ev) = (eps(0), eps(1));
|
||||
let same = |a: &[f64; 2], b: &[f64; 2]| (a[0] - b[0]).abs() <= eu && (a[1] - b[1]).abs() <= ev;
|
||||
let (mut uv, mut pts) = (Vec::with_capacity(l.uv.len()), Vec::with_capacity(l.pts.len()));
|
||||
for (q, p) in l.uv.iter().zip(&l.pts) {
|
||||
if uv.last().is_some_and(|last| same(last, q)) {
|
||||
continue;
|
||||
}
|
||||
uv.push(*q);
|
||||
pts.push(*p);
|
||||
}
|
||||
// The loop is cyclic: its last node may land on its first.
|
||||
while uv.len() > 2 && uv.last().zip(uv.first()).is_some_and(|(a, b)| same(a, b)) {
|
||||
uv.pop();
|
||||
pts.pop();
|
||||
}
|
||||
l.uv = uv;
|
||||
l.pts = pts;
|
||||
}
|
||||
|
||||
/// Bring every loop of a face into the same period as its outer boundary. A loop is projected
|
||||
/// on its own, and a periodic surface repeats, so a hole can come out a whole turn away from the
|
||||
/// boundary it lies in — where it no longer cuts the face, and the face is drawn over its hole.
|
||||
fn align_periods(loops: &mut [seam::Loop], dom: &uv::Domain) {
|
||||
for d in 0..2 {
|
||||
let Some(per) = dom.period.get(d).copied().flatten().filter(|p| *p > 0.0) else { continue };
|
||||
let middle = |l: &seam::Loop| {
|
||||
let (lo, hi) = l.uv.iter().filter_map(|q| q.get(d)).fold((f64::MAX, f64::MIN), |(a, b), x| (a.min(*x), b.max(*x)));
|
||||
0.5 * (lo + hi)
|
||||
};
|
||||
let Some(outer) = loops.first().map(middle) else { continue };
|
||||
if !outer.is_finite() {
|
||||
continue;
|
||||
}
|
||||
for l in loops.iter_mut().skip(1) {
|
||||
let shift = ((outer - middle(l)) / per).round() * per;
|
||||
if shift != 0.0 && shift.is_finite() {
|
||||
for q in l.uv.iter_mut() {
|
||||
if let Some(x) = q.get_mut(d) {
|
||||
*x += shift;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Nodes a face's boundary needs before it can hold a triangle.
|
||||
const MIN_BOUNDARY_NODES: usize = 3;
|
||||
|
||||
/// A face whose whole boundary comes to fewer than three nodes has no triangle to stand on: it is
|
||||
/// a sliver no wider than the tolerance (a chord with the shallow arc beside it, which the edge
|
||||
/// sampler rightly leaves as one segment each). Its edges are sampled again, halved until the
|
||||
/// boundary has nodes enough, so the sliver is drawn instead of failing the whole body's plan.
|
||||
/// The nodes stay shared, so the faces beside it follow the same points.
|
||||
fn divide_collapsed_faces(faces: &[&mt::Face], jobs: &[(&mt::Edge, f64)], rules: &size::Rules, nodes: &mut HashMap<mt::EdgeID, edges::Nodes>) {
|
||||
/// Halvings before giving up (a boundary of two edges reaches four nodes in one).
|
||||
const MAX_HALVINGS: usize = 8;
|
||||
let h_of: HashMap<mt::EdgeID, (&mt::Edge, f64)> = jobs.iter().map(|(e, h)| (e.id(), (*e, *h))).collect();
|
||||
for f in faces {
|
||||
let ids: Vec<mt::EdgeID> = f.absolute_boundaries().iter().flat_map(|w| w.edge_iter().map(|e| e.id())).collect();
|
||||
let count = |nodes: &HashMap<mt::EdgeID, edges::Nodes>| {
|
||||
ids.iter().map(|id| nodes.get(id).map_or(0, |n| n.points.len().saturating_sub(1))).sum::<usize>()
|
||||
};
|
||||
if count(nodes) >= MIN_BOUNDARY_NODES {
|
||||
continue;
|
||||
}
|
||||
for _ in 0..MAX_HALVINGS {
|
||||
for id in &ids {
|
||||
let Some((e, h)) = h_of.get(id) else { continue };
|
||||
let longest =
|
||||
nodes.get(id).map_or(0.0, |n| n.points.windows(2).filter_map(|w| Some(w.first()?.dist(*w.get(1)?))).fold(0.0, f64::max));
|
||||
let finer = if longest > 0.0 { (longest * 0.5).min(*h) } else { *h };
|
||||
nodes.insert(*id, edges::sample_edge(e, rules, finer));
|
||||
}
|
||||
if count(nodes) >= MIN_BOUNDARY_NODES {
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// How far a boundary node may lie off its surface and still be projected onto it (see the note
|
||||
/// in [`plan_face`]).
|
||||
fn off_surface(rules: &size::Rules) -> f64 {
|
||||
(rules.extent * ON_SURFACE).max(rules.boundary_tol.unwrap_or(0.0)).max(1e-6)
|
||||
}
|
||||
|
||||
/// Whether a face of `solid` has a loop that runs along a seam edge between two distinct
|
||||
/// vertices and straight back along it, the two uses one after the other (they meet at a pole):
|
||||
/// a whole sphere as one face bounded by its meridian alone, as SMLib writes it, or a sphere
|
||||
/// face whose meridian is split where a trim crosses it. Truck's mesher projects both uses onto
|
||||
/// one parameter line and meshes nothing; the display plan opens the loop across the period
|
||||
/// ([`open_pole_to_pole`], and [`open_at_pole`] for a seam from one pole to a rim).
|
||||
pub(crate) fn has_pole_to_pole_seam(solid: &Solid) -> bool {
|
||||
solid.face_iter().any(|f| {
|
||||
f.boundaries().iter().any(|w| {
|
||||
let edges: Vec<&mt::Edge> = w.edge_iter().collect();
|
||||
let n = edges.len();
|
||||
n >= 2
|
||||
&& (0..n).any(|k| match (edges.get(k), edges.get((k + 1) % n)) {
|
||||
(Some(a), Some(b)) => a.id() == b.id() && a.orientation() != b.orientation() && a.front() != a.back(),
|
||||
_ => false,
|
||||
})
|
||||
})
|
||||
})
|
||||
Ok(())
|
||||
}
|
||||
|
||||
/// Twice the signed area of a parameter polygon.
|
||||
@@ -636,9 +411,9 @@ fn holes_beside_outer(dom: &uv::Domain, loops: &mut [seam::Loop]) {
|
||||
/// `size`: the face's target side length (mm) and the surface's mean speed along u and v.
|
||||
fn plan_face(
|
||||
fi: usize,
|
||||
kind: &str,
|
||||
s: &EvalSurface,
|
||||
loops3: &[Vec<Vec3>],
|
||||
known: &[Vec<Option<[f64; 2]>>],
|
||||
rules: &size::Rules,
|
||||
size: (f64, [f64; 2]),
|
||||
stats: &mut PlanStats,
|
||||
@@ -646,26 +421,17 @@ fn plan_face(
|
||||
let tol = rules.tol;
|
||||
let (h, speed) = size;
|
||||
let dom = uv::domain(s);
|
||||
// How far a boundary node may lie off its surface and still be projected onto it: a property
|
||||
// of the body, not of the mesh asked for, so a finer chord tolerance never refuses a boundary
|
||||
// a coarser one accepted, and zooming in never drops a face (#1721). A hundredth of the body,
|
||||
// or more where the file it was read from declares more; a file's declared accuracy never
|
||||
// lowers it, since writers' edges stray farther from their surfaces than they declare.
|
||||
let accept = off_surface(rules);
|
||||
let fail = |why: &str| unsupported(format!("face {fi}: {why}"));
|
||||
let accept = (tol * 10.0).max(1e-6);
|
||||
let fail = |why: &str| unsupported(format!("face {fi} ({kind}): {why}"));
|
||||
// Project every loop, unwrapping periodic parameters along it.
|
||||
let mut projected: Vec<seam::Loop> = Vec::with_capacity(loops3.len());
|
||||
for (li, pts) in loops3.iter().enumerate() {
|
||||
let mut luv: Vec<[f64; 2]> = Vec::with_capacity(pts.len());
|
||||
let mut prev: Option<[f64; 2]> = None;
|
||||
for (k, p) in pts.iter().enumerate() {
|
||||
// Known parameters (an intersection edge's), else a projection seeded by the
|
||||
// previous node.
|
||||
let given = known.get(li).and_then(|l| l.get(k)).copied().flatten();
|
||||
let mut q = match given {
|
||||
Some(q) => q,
|
||||
None => uv::project(s, &dom, *p, prev, accept).ok_or_else(|| fail("boundary does not lie on the surface"))?,
|
||||
};
|
||||
// A projection seeded by the previous node.
|
||||
let mut q = uv::project(s, &dom, *p, prev, accept)
|
||||
.ok_or_else(|| fail(&format!("boundary node {k} of loop {li} at {p:?} does not lie on the surface")))?;
|
||||
if let Some(pr) = prev {
|
||||
for d in 0..2 {
|
||||
q[d] = uv::unwrap(q[d], pr[d], dom.period[d]);
|
||||
|
||||
@@ -384,118 +384,3 @@ fn seam_nodes(s: &EvalSurface, swap: bool, x: f64, y0: f64, y1: f64, rules: &Rul
|
||||
let inner = ys.get(1..ys.len() - 1)?;
|
||||
Some(inner.iter().map(|&y| (y, s.eval(uv(y)[0], uv(y)[1]).0)).collect())
|
||||
}
|
||||
|
||||
/// Whether a face of `solid` has a boundary loop turning more than once around the axis its
|
||||
/// surface turns about (a helical trim): a surface of revolution's own axis, or the axis a
|
||||
/// spline's iso-line winds around (the flank of a screw sweep, as SMLib writes a modelled
|
||||
/// thread's groove, whose bore or shaft it splits into one face per turn). Truck's mesher cannot
|
||||
/// mesh such a face (its triangles cut across the solid, or its triangulation fails), so the
|
||||
/// reference mesh is made from the display plan instead.
|
||||
pub(crate) fn turns_more_than_once(solid: &crate::body::Solid) -> bool {
|
||||
use truck_meshalgo::prelude::*;
|
||||
use truck_modeling as mt;
|
||||
let tau = std::f64::consts::TAU;
|
||||
for face in solid.face_iter() {
|
||||
let Some((o, k)) = winding_axis(&face.surface()) else { continue };
|
||||
let e1 = k.cross(if k.x.abs() < 0.9 { Vec3::X } else { Vec3::Y }).normalized().unwrap_or(Vec3::X);
|
||||
let e2 = k.cross(e1);
|
||||
let angle = |p: mt::Point3| {
|
||||
let d = crate::body::from_p3(p) - o;
|
||||
d.dot(e2).atan2(d.dot(e1))
|
||||
};
|
||||
for wire in face.boundaries() {
|
||||
let (mut total, mut lo, mut hi) = (0.0f64, 0.0f64, 0.0f64);
|
||||
let mut prev: Option<f64> = None;
|
||||
for e in wire.edge_iter() {
|
||||
let c = e.oriented_curve();
|
||||
let (t0, t1) = c.range_tuple();
|
||||
// Samples dense enough that no step turns half way round.
|
||||
let mut n = 16usize;
|
||||
let steps = loop {
|
||||
let xs: Vec<f64> = (0..=n).map(|i| angle(c.subs(t0 + (t1 - t0) * i as f64 / n as f64))).collect();
|
||||
let ok = xs.windows(2).all(|w| half_turn(w[1] - w[0]).abs() < 1.0);
|
||||
if ok || n >= 4096 {
|
||||
break xs;
|
||||
}
|
||||
n *= 4;
|
||||
};
|
||||
for x in steps {
|
||||
if let Some(p) = prev {
|
||||
total += half_turn(x - p);
|
||||
lo = lo.min(total);
|
||||
hi = hi.max(total);
|
||||
}
|
||||
prev = Some(x);
|
||||
}
|
||||
}
|
||||
if hi - lo > tau * 1.01 {
|
||||
return true;
|
||||
}
|
||||
}
|
||||
}
|
||||
false
|
||||
}
|
||||
|
||||
/// An angle step taken the short way round, in (−π, π].
|
||||
fn half_turn(x: f64) -> f64 {
|
||||
let pi = std::f64::consts::PI;
|
||||
(x + pi).rem_euclid(std::f64::consts::TAU) - pi
|
||||
}
|
||||
|
||||
/// Samples per knot span of a spline's iso-line in [`winding_axis`], and their bounds.
|
||||
const ISO_SPAN_SAMPLES: usize = 4;
|
||||
const ISO_MIN_SAMPLES: usize = 64;
|
||||
const ISO_MAX_SAMPLES: usize = 8192;
|
||||
|
||||
/// The axis (a point on it, its unit direction) a face's surface turns about, if it turns at
|
||||
/// all: a surface of revolution's own axis; for a spline, the axis an iso-line through the
|
||||
/// middle of the other parameter winds more than once around (a screw sweep's flank, a coil's
|
||||
/// wall). That axis is the mean binormal of the iso-line (each pair of successive tangents
|
||||
/// turns about it), through the mean of its samples; `None` for a spline neither iso-line of
|
||||
/// which turns more than once about it.
|
||||
fn winding_axis(surface: &truck_modeling::Surface) -> Option<(Vec3, Vec3)> {
|
||||
use truck_meshalgo::prelude::*;
|
||||
use truck_modeling as mt;
|
||||
if let mt::Surface::RevolutedCurve(proc) = surface {
|
||||
let rev = proc.entity();
|
||||
let m = *proc.transform();
|
||||
let o = crate::body::from_p3(m.transform_point(rev.origin()));
|
||||
let a = m.transform_vector(rev.axis());
|
||||
return Vec3::new(a.x, a.y, a.z).normalized().map(|k| (o, k));
|
||||
}
|
||||
if matches!(surface, mt::Surface::Plane(_)) {
|
||||
return None;
|
||||
}
|
||||
let s = super::surface::eval_surface(surface)?;
|
||||
let EvalSurface::Nurbs { knots, .. } = &s else { return None };
|
||||
let dom = super::uv::domain(&s);
|
||||
let tau = std::f64::consts::TAU;
|
||||
(0..2).find_map(|d| {
|
||||
let (a, b) = dom.range[d]?;
|
||||
let (c0, c1) = dom.range[1 - d]?;
|
||||
let across = 0.5 * (c0 + c1);
|
||||
let n = knots.get(d).map_or(0, |k| k.len()).saturating_mul(ISO_SPAN_SAMPLES).clamp(ISO_MIN_SAMPLES, ISO_MAX_SAMPLES);
|
||||
let pts: Vec<Vec3> = (0..=n)
|
||||
.filter_map(|i| {
|
||||
let t = a + (b - a) * i as f64 / n as f64;
|
||||
s.derivs(if d == 0 { t } else { across }, if d == 0 { across } else { t }).map(|(p, _, _)| p)
|
||||
})
|
||||
.collect();
|
||||
if pts.len() != n + 1 {
|
||||
return None;
|
||||
}
|
||||
let tangents: Vec<Vec3> = pts.windows(2).filter_map(|w| (w[1] - w[0]).normalized()).collect();
|
||||
let binormal = tangents.windows(2).fold(Vec3::ZERO, |acc, w| acc + w[0].cross(w[1]));
|
||||
let k = binormal.normalized()?;
|
||||
let o = pts.iter().fold(Vec3::ZERO, |acc, p| acc + *p) * (1.0 / pts.len() as f64);
|
||||
let e1 = k.cross(if k.x.abs() < 0.9 { Vec3::X } else { Vec3::Y }).normalized()?;
|
||||
let e2 = k.cross(e1);
|
||||
let angles: Vec<f64> = pts.iter().map(|p| (*p - o).dot(e2).atan2((*p - o).dot(e1))).collect();
|
||||
// Too coarse to follow the turns (a step of a radian or more): not taken as winding.
|
||||
let steps: Vec<f64> = angles.windows(2).map(|w| half_turn(w[1] - w[0])).collect();
|
||||
if steps.iter().any(|x| x.abs() >= 1.0) {
|
||||
return None;
|
||||
}
|
||||
(steps.iter().sum::<f64>().abs() > tau * 1.01).then_some((o, k))
|
||||
})
|
||||
}
|
||||
|
||||
@@ -26,6 +26,13 @@ const MIN_LEN_FRACTION: f64 = 1.0 / 4.0;
|
||||
/// Curvature samples per parameter direction.
|
||||
const SAMPLES: usize = 9;
|
||||
|
||||
/// Most segments a face's boundary may get from the facet angle: past it, the face's curved
|
||||
/// edges are sampled again at twice the angle (a plate with a thousand small holes is not worth
|
||||
/// a hundred times the work to draw each hole round).
|
||||
pub(crate) const MAX_ANGLE_BOUNDARY: usize = 4096;
|
||||
/// How often the angle is doubled for a crowded face before only the chord tolerance is left.
|
||||
pub(crate) const MAX_RELAX: u8 = 4;
|
||||
|
||||
/// The most one facet may turn (radians): an edge segment's tangent or a cell's surface normal.
|
||||
/// The angle divides only where a facet would sag more than `floor` (mm), so a radius far below
|
||||
/// the tolerance keeps its chord facets.
|
||||
@@ -41,12 +48,11 @@ impl Turn {
|
||||
sag > self.floor && angle > self.angle
|
||||
}
|
||||
|
||||
/// The angle doubled `level` times, as the reference mesher relaxes it for a crowded face;
|
||||
/// from its last level, or past half a turn, no angle at all (the tolerance alone).
|
||||
/// The angle doubled `level` times, for a crowded face (see [`MAX_ANGLE_BOUNDARY`]); from
|
||||
/// level [`MAX_RELAX`], or past half a turn, no angle at all (the tolerance alone).
|
||||
pub(crate) fn relaxed(&self, level: u8) -> Option<Turn> {
|
||||
use truck_meshalgo::tessellation::Deflection;
|
||||
let angle = self.angle * f64::from(1u32 << u32::from(level.min(Deflection::MAX_RELAX)));
|
||||
(level == 0 || (level < Deflection::MAX_RELAX && angle < std::f64::consts::PI)).then_some(Turn { angle, ..*self })
|
||||
let angle = self.angle * f64::from(1u32 << u32::from(level.min(MAX_RELAX)));
|
||||
(level == 0 || (level < MAX_RELAX && angle < std::f64::consts::PI)).then_some(Turn { angle, ..*self })
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@@ -1,282 +0,0 @@
|
||||
//! Points on a surface–surface intersection edge, on our own surface evaluators.
|
||||
//!
|
||||
//! truck stores an intersection edge as a leader curve (a polyline or B-spline near the
|
||||
//! intersection) and its two surfaces, and finds each point by projecting the leader point onto
|
||||
//! each surface in turn from no starting guess: a nearest-point search per surface per point.
|
||||
//! Here the intersection is solved exactly at a set of leader parameters, and a cubic B-spline
|
||||
//! is fitted through those points together with both surfaces' parameters there (a smooth
|
||||
//! leader with its curves on both surfaces). The B-spline runs over the chord length of the
|
||||
//! solved points, never over the curve's own parameter: truck's parameter can change speed
|
||||
//! abruptly (a polyline leader with short and long spans), and a cubic interpolated over it
|
||||
//! rings there, stepping back along the curve (Lee, "Choosing nodes in parametric curve
|
||||
//! interpolation", CAD 21, 1989). The curve's parameter reaches the chord length through the
|
||||
//! piecewise-linear map between the seeds, which is monotone. Each point is then refined from
|
||||
//! the B-spline by Newton's method on the four unknowns (u1, v1, u2, v2): S1(u1, v1) = S2(u2, v2), held on the
|
||||
//! plane through the B-spline point normal to the B-spline (Barnhill & Kersey, "A marching
|
||||
//! method for parametric surface/surface intersection", CAGD 7, 1990; Patrikalakis & Maekawa,
|
||||
//! "Shape Interrogation for CAD/CAM", 2002, ch. 5). Every point is independent once seeded. Where the surfaces are close to
|
||||
//! tangent along the edge (the system becomes singular) or Newton does not converge, the caller
|
||||
//! uses truck's evaluation.
|
||||
|
||||
use solidcraft_geom::Vec3;
|
||||
use solidcraft_geom::bspline::CurveN;
|
||||
use solidcraft_geom::tess::EvalSurface;
|
||||
use truck_meshalgo::prelude::*;
|
||||
use truck_modeling as mt;
|
||||
|
||||
use super::uv::{self, Domain};
|
||||
use crate::body::from_p3;
|
||||
|
||||
/// Leader parameters at which both surfaces' parameters are found up front.
|
||||
const SEEDS: usize = 32;
|
||||
/// Newton iterations per point (two or three suffice from a seed).
|
||||
const MAX_ITER: usize = 8;
|
||||
/// Below this sine of the angle between the surface normals the intersection is treated as
|
||||
/// near-tangent: points there are left to truck.
|
||||
const MIN_SIN: f64 = 0.05;
|
||||
|
||||
/// An intersection edge ready for point evaluation.
|
||||
pub(crate) struct Intersection {
|
||||
surfaces: [EvalSurface; 2],
|
||||
domains: [Domain; 2],
|
||||
/// The smooth leader over the chord length of the seeds: the intersection point (x, y, z)
|
||||
/// and both surfaces' parameters (u1, v1, u2, v2), periodic directions unwrapped along it.
|
||||
leader: CurveN<7>,
|
||||
/// The curve's parameter at each seed and the chord length there, both increasing.
|
||||
seeds: Vec<(f64, f64)>,
|
||||
/// Convergence tolerance (mm).
|
||||
eps: f64,
|
||||
}
|
||||
|
||||
impl Intersection {
|
||||
/// `None` unless `curve` is an intersection curve whose surfaces have evaluators and whose
|
||||
/// leader projects onto both.
|
||||
pub(crate) fn new(curve: &mt::Curve) -> Option<Self> {
|
||||
let mt::Curve::IntersectionCurve(ic) = curve else { return None };
|
||||
let surfaces = [super::surface::eval_surface(ic.surface0())?, super::surface::eval_surface(ic.surface1())?];
|
||||
let domains = [uv::domain(&surfaces[0]), uv::domain(&surfaces[1])];
|
||||
let (t0, t1) = curve.range_tuple();
|
||||
if !(t0.is_finite() && t1.is_finite() && t1 > t0) {
|
||||
return None;
|
||||
}
|
||||
let leader = |t: f64| match ic.leader() {
|
||||
mt::Leader::Polyline(l) => (from_p3(l.subs(t)), from_v3(l.der(t))),
|
||||
mt::Leader::BSpline(l) => (from_p3(l.subs(t)), from_v3(l.der(t))),
|
||||
};
|
||||
// The leader lies within the curve's tolerance of both surfaces.
|
||||
let accept = (ic.tolerance() * 100.0).max(1e-6);
|
||||
let scale = (0..=SEEDS).map(|k| leader(t0 + (t1 - t0) * k as f64 / SEEDS as f64).0.len()).fold(1.0, f64::max);
|
||||
let mut it = Intersection {
|
||||
surfaces,
|
||||
domains,
|
||||
leader: CurveN { degree: 1, knots: Vec::new(), ctrl: Vec::new() },
|
||||
seeds: Vec::new(),
|
||||
eps: scale * 1e-10,
|
||||
};
|
||||
let (mut params, mut points): (Vec<f64>, Vec<[f64; 7]>) = (Vec::with_capacity(SEEDS + 1), Vec::with_capacity(SEEDS + 1));
|
||||
let mut prev: Option<[f64; 4]> = None;
|
||||
for k in 0..=SEEDS {
|
||||
let t = t0 + (t1 - t0) * k as f64 / SEEDS as f64;
|
||||
let (p, d) = leader(t);
|
||||
let mut q = [0.0; 4];
|
||||
for (i, (s, dom)) in it.surfaces.iter().zip(&it.domains).enumerate() {
|
||||
let hint = prev.map(|w| [w[2 * i], w[2 * i + 1]]);
|
||||
let mut uv = uv::project(s, dom, p, hint, accept)?;
|
||||
if let Some(h) = hint {
|
||||
for d in 0..2 {
|
||||
uv[d] = uv::unwrap(uv[d], h[d], dom.period[d]);
|
||||
}
|
||||
}
|
||||
q[2 * i] = uv[0];
|
||||
q[2 * i + 1] = uv[1];
|
||||
}
|
||||
// The fit takes only points exactly on the intersection: a leader point Newton
|
||||
// cannot place (near-tangent, or the solution leaves a surface's domain where a
|
||||
// polyline leader overshoots a vertex) is off the curve and would bend the fit back
|
||||
// along it. The ends stay either way, so the fit spans the whole curve.
|
||||
let solved = d.normalized().and_then(|dir| it.newton(p, dir, q));
|
||||
prev = Some(solved.map_or(q, |(_, q)| q));
|
||||
let (p, q) = match solved {
|
||||
Some(x) => x,
|
||||
None if k == 0 || k == SEEDS => (p, q),
|
||||
None => continue,
|
||||
};
|
||||
// A seed that lands on the previous one adds no chord and stays out of the fit
|
||||
// (whose parameters must increase).
|
||||
let s = match it.seeds.last().zip(points.last()) {
|
||||
Some((&(_, s), &[x, y, z, ..])) => {
|
||||
let chord = p.dist(Vec3::new(x, y, z));
|
||||
if !(chord > it.eps) {
|
||||
continue;
|
||||
}
|
||||
s + chord
|
||||
}
|
||||
None => 0.0,
|
||||
};
|
||||
it.seeds.push((t, s));
|
||||
params.push(s);
|
||||
points.push([p.x, p.y, p.z, q[0], q[1], q[2], q[3]]);
|
||||
}
|
||||
let degree = points.len().checked_sub(1)?.min(3);
|
||||
it.leader = CurveN::interpolate(¶ms, &points, degree).ok()?;
|
||||
Some(it)
|
||||
}
|
||||
|
||||
/// The point at curve parameter `t` and both surfaces' parameters there, or `None`
|
||||
/// (near-tangent, or Newton failed).
|
||||
pub(crate) fn eval(&self, t: f64) -> Option<(Vec3, [f64; 4])> {
|
||||
let (x, dx) = self.leader.eval(self.chord_at(t)?)?;
|
||||
let dir = Vec3::new(dx[0], dx[1], dx[2]).normalized()?;
|
||||
self.newton(Vec3::new(x[0], x[1], x[2]), dir, [x[3], x[4], x[5], x[6]])
|
||||
}
|
||||
|
||||
/// The chord length at curve parameter `t`: linear between the seeds, clamped to the ends.
|
||||
fn chord_at(&self, t: f64) -> Option<f64> {
|
||||
let (&(ta, sa), &(tb, sb)) = (self.seeds.first()?, self.seeds.last()?);
|
||||
if !(t > ta) {
|
||||
return Some(sa);
|
||||
}
|
||||
if !(t < tb) {
|
||||
return Some(sb);
|
||||
}
|
||||
let k = self.seeds.partition_point(|&(tk, _)| tk <= t);
|
||||
let (Some(&(t0, s0)), Some(&(t1, s1))) = (<[_]>::get(&self.seeds, k.checked_sub(1)?), <[_]>::get(&self.seeds, k)) else { return None };
|
||||
Some(s0 + (s1 - s0) * (t - t0) / (t1 - t0))
|
||||
}
|
||||
|
||||
/// The two surfaces, in the order of the parameters [`Intersection::eval`] gives.
|
||||
pub(crate) fn surfaces(&self) -> &[EvalSurface; 2] {
|
||||
&self.surfaces
|
||||
}
|
||||
|
||||
/// Newton from the parameters `x` on S1 = S2, held on the plane through `l` normal to `dir`:
|
||||
/// the point and its parameters, or `None` (near-tangent, singular, or not converged).
|
||||
fn newton(&self, l: Vec3, dir: Vec3, mut x: [f64; 4]) -> Option<(Vec3, [f64; 4])> {
|
||||
for _ in 0..MAX_ITER {
|
||||
let (Some((p1, u1, v1)), Some((p2, u2, v2))) = (self.surfaces[0].derivs(x[0], x[1]), self.surfaces[1].derivs(x[2], x[3])) else {
|
||||
return None;
|
||||
};
|
||||
let gap = p1 - p2;
|
||||
let plane = dir.dot(p1 - l);
|
||||
let n = u1.cross(v1).normalized()?.cross(u2.cross(v2).normalized()?);
|
||||
if n.len() < MIN_SIN {
|
||||
return None;
|
||||
}
|
||||
if gap.len() <= self.eps && plane.abs() <= self.eps {
|
||||
return Some(((p1 + p2) * 0.5, x));
|
||||
}
|
||||
// J · dx = −F, F = (S1 − S2, dir · (S1 − L)).
|
||||
let j = [[u1.x, v1.x, -u2.x, -v2.x], [u1.y, v1.y, -u2.y, -v2.y], [u1.z, v1.z, -u2.z, -v2.z], [dir.dot(u1), dir.dot(v1), 0.0, 0.0]];
|
||||
let dx = solve4(j, [-gap.x, -gap.y, -gap.z, -plane])?;
|
||||
for (xi, d) in x.iter_mut().zip(dx) {
|
||||
*xi += d;
|
||||
}
|
||||
for i in 0..2 {
|
||||
let [u, v] = self.domains[i].clamp([x[2 * i], x[2 * i + 1]]);
|
||||
x[2 * i] = u;
|
||||
x[2 * i + 1] = v;
|
||||
}
|
||||
}
|
||||
None
|
||||
}
|
||||
}
|
||||
|
||||
fn from_v3(v: mt::Vector3) -> Vec3 {
|
||||
Vec3::new(v.x, v.y, v.z)
|
||||
}
|
||||
|
||||
/// `a · x = b` for a 4 × 4 system by Gaussian elimination with partial pivoting; `None` when
|
||||
/// singular.
|
||||
fn solve4(mut a: [[f64; 4]; 4], mut b: [f64; 4]) -> Option<[f64; 4]> {
|
||||
let scale = a.iter().flatten().fold(0.0f64, |m, x| m.max(x.abs()));
|
||||
if !(scale > 0.0) {
|
||||
return None;
|
||||
}
|
||||
for c in 0..4 {
|
||||
let p = (c..4).max_by(|&i, &j| a[i][c].abs().total_cmp(&a[j][c].abs()))?;
|
||||
if !(a[p][c].abs() > 1e-12 * scale) {
|
||||
return None;
|
||||
}
|
||||
a.swap(c, p);
|
||||
b.swap(c, p);
|
||||
let pivot = a[c];
|
||||
for r in c + 1..4 {
|
||||
let f = a[r][c] / pivot[c];
|
||||
for (x, p) in a[r].iter_mut().zip(pivot).skip(c) {
|
||||
*x -= f * p;
|
||||
}
|
||||
b[r] -= f * b[c];
|
||||
}
|
||||
}
|
||||
let mut x = [0.0; 4];
|
||||
for c in (0..4).rev() {
|
||||
let s: f64 = (c + 1..4).map(|k| a[c][k] * x[k]).sum();
|
||||
x[c] = (b[c] - s) / a[c][c];
|
||||
}
|
||||
x.iter().all(|v| v.is_finite()).then_some(x)
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
#[test]
|
||||
fn solves_a_four_by_four_system() {
|
||||
let a = [[2.0, 1.0, 0.0, 0.0], [1.0, 3.0, 1.0, 0.0], [0.0, 1.0, 4.0, 1.0], [0.0, 0.0, 1.0, 5.0]];
|
||||
let x = [1.0, -2.0, 0.5, 3.0];
|
||||
let b = [0, 1, 2, 3].map(|r| (0..4).map(|c| a[r][c] * x[c]).sum::<f64>());
|
||||
let Some(got) = solve4(a, b) else { panic!("solves") };
|
||||
for (g, w) in got.iter().zip(x) {
|
||||
assert!((g - w).abs() < 1e-12);
|
||||
}
|
||||
assert!(solve4([[1.0; 4]; 4], [1.0; 4]).is_none(), "singular");
|
||||
}
|
||||
|
||||
/// Points on the intersection edges of a boolean (a cross hole: two cylinders) lie on both
|
||||
/// surfaces, on truck's curve (the parameterisation follows our smooth leader, so the same
|
||||
/// parameter may sit at a slightly different place along it), and advance along it.
|
||||
#[test]
|
||||
fn intersection_points_lie_on_both_surfaces_and_on_truck_s_curve() {
|
||||
let Ok(pin) = crate::cylinder(Vec3::new(5.0, 5.0, -1.0), Vec3::Z, 2.0, 6.0) else { panic!("pin") };
|
||||
let Ok(cross) = crate::cylinder(Vec3::new(5.0, -1.0, 2.0), Vec3::Y, 1.0, 12.0) else { panic!("cross") };
|
||||
let Ok(Some(body)) = crate::boolean(&pin, &cross, crate::BoolOp::Cut) else { panic!("cut") };
|
||||
let mut checked = 0;
|
||||
for e in body.solid.edge_iter() {
|
||||
let curve = e.curve();
|
||||
let Some(ix) = Intersection::new(&curve) else { continue };
|
||||
let (t0, t1) = curve.range_tuple();
|
||||
let dense: Vec<Vec3> = (0..=1000).map(|k| from_p3(curve.subs(t0 + (t1 - t0) * f64::from(k) / 1000.0))).collect();
|
||||
let near = |p: Vec3| dense.windows(2).map(|w| p.dist_to_segment(w[0], w[1])).fold(f64::MAX, f64::min);
|
||||
// Arc length along truck's curve to the point nearest `p`.
|
||||
let along = |p: Vec3| {
|
||||
let (mut run, mut best) = (0.0, (f64::MAX, 0.0));
|
||||
for w in dense.windows(2) {
|
||||
let (a, b) = (w[0], w[1]);
|
||||
let len = a.dist(b);
|
||||
let f = if len > 0.0 { ((p - a).dot(b - a) / (len * len)).clamp(0.0, 1.0) } else { 0.0 };
|
||||
let d = p.dist(a + (b - a) * f);
|
||||
if d < best.0 {
|
||||
best = (d, run + f * len);
|
||||
}
|
||||
run += len;
|
||||
}
|
||||
best.1
|
||||
};
|
||||
let mut last: Option<Vec3> = None;
|
||||
for k in 0..=50 {
|
||||
let t = t0 + (t1 - t0) * f64::from(k) / 50.0;
|
||||
let Some((p, _)) = ix.eval(t) else { continue };
|
||||
for (s, dom) in ix.surfaces.iter().zip(&ix.domains) {
|
||||
assert!(uv::project(s, dom, p, None, 1e-8).is_some(), "t {t}: {p:?} off a surface");
|
||||
}
|
||||
assert!(near(p) < 1e-4, "t {t}: {p:?} is {} from truck's curve", near(p));
|
||||
if let Some(q) = last {
|
||||
assert!(along(p) > along(q), "t {t}: the points advance along the curve");
|
||||
}
|
||||
last = Some(p);
|
||||
checked += 1;
|
||||
}
|
||||
}
|
||||
assert!(checked > 100, "{checked} points checked");
|
||||
}
|
||||
}
|
||||
+103
-110
@@ -1,119 +1,112 @@
|
||||
//! Truck surfaces in the evaluator's form ([`EvalSurface`]), which the CPU and GPU evaluators
|
||||
//! share. Surfaces the evaluators cannot express (a revolved intersection curve) or that are
|
||||
//! too large for them give `None`; their faces are evaluated here instead.
|
||||
//! Surface records ([`brep::Surface`]) in the evaluator's form ([`EvalSurface`]), which the
|
||||
//! planner and the CPU and GPU evaluators share. Planes, splines and revolutions keep their
|
||||
//! own parameters ([`SurfaceEval`] holds them in that form already). A cylinder, cone, sphere
|
||||
//! or torus is the revolution it is, exactly: its generator (a line, or a circle as an exact
|
||||
//! rational arc) turned about its axis, the angle as the revolution's parameter. The planner
|
||||
//! finds a boundary's parameters by projecting its nodes, so the generator's own parameter
|
||||
//! need not be the record's height or latitude.
|
||||
|
||||
use solidcraft_geom::Vec3;
|
||||
use solidcraft_geom::bspline::MAX_SPLINE_DEGREE;
|
||||
use solidcraft_geom::tess::{EvalCurve, EvalSurface, MAX_EVAL_KNOTS, MAX_EVAL_POINTS};
|
||||
use truck_meshalgo::prelude::*;
|
||||
use truck_modeling as mt;
|
||||
use solidcraft_geom::brep::eval::SurfaceEval;
|
||||
use solidcraft_geom::brep::{self, Frame};
|
||||
use solidcraft_geom::tess::{EvalCurve, EvalSurface};
|
||||
|
||||
use crate::body::from_p3;
|
||||
use crate::{KernelError, Result};
|
||||
|
||||
/// `surface` in the evaluator's form, parameterised exactly as truck parameterises it.
|
||||
pub(crate) fn eval_surface(surface: &mt::Surface) -> Option<EvalSurface> {
|
||||
let s = match surface {
|
||||
mt::Surface::Plane(p) => {
|
||||
EvalSurface::Plane { origin: from_p3(p.subs(0.0, 0.0)), du: from_v3(p.uder(0.0, 0.0)), dv: from_v3(p.vder(0.0, 0.0)) }
|
||||
}
|
||||
mt::Surface::BSplineSurface(b) => {
|
||||
let ctrl = b.control_points().iter().map(|row| row.iter().map(|p| [p.x, p.y, p.z, 1.0]).collect()).collect();
|
||||
nurbs(b.degrees(), b.knot_vecs(), ctrl)?
|
||||
}
|
||||
mt::Surface::NurbsSurface(n) => {
|
||||
let ctrl = n.control_points().iter().map(|row| row.iter().map(|h| [h.x, h.y, h.z, h.w]).collect()).collect();
|
||||
nurbs((n.udegree(), n.vdegree()), n.knot_vecs(), ctrl)?
|
||||
}
|
||||
mt::Surface::RevolutedCurve(proc) => {
|
||||
let rev = proc.entity();
|
||||
let profile = eval_curve(rev.entity_curve())?;
|
||||
let axis = from_v3(rev.axis()).normalized()?;
|
||||
let m = proc.transform();
|
||||
// cgmath is column-major: m[column][row].
|
||||
let transform = [0, 1, 2].map(|r| [m[0][r], m[1][r], m[2][r], m[3][r]]);
|
||||
if transform.iter().flatten().any(|x| !x.is_finite()) {
|
||||
return None;
|
||||
}
|
||||
EvalSurface::Revolved { profile, origin: from_p3(rev.origin()), axis, transform, swap: !proc.orientation() }
|
||||
}
|
||||
fn v3(p: [f64; 3]) -> Vec3 {
|
||||
Vec3::new(p[0], p[1], p[2])
|
||||
}
|
||||
|
||||
/// No transform: the revolution as it is.
|
||||
const IDENTITY: [[f64; 4]; 3] = [[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0]];
|
||||
|
||||
/// `surface` in the evaluator's form. `domain` is the face's parameter box in the record's
|
||||
/// parameters, `[[u0, u1], [v0, v1]]`: a cylinder's or cone's generator spans the face's
|
||||
/// heights (an unbounded line has no parameter box to project into).
|
||||
pub(crate) fn eval_surface(surface: &brep::Surface, domain: [[f64; 2]; 2]) -> Result<EvalSurface> {
|
||||
use std::f64::consts::{FRAC_PI_2, TAU};
|
||||
let invalid = |what: &str| KernelError::Invalid(format!("tessellation: {what}"));
|
||||
let axes = |f: &Frame| {
|
||||
let (o, x, y) = (v3(f.origin), v3(f.x_axis), v3(f.y_axis));
|
||||
(o, x, x.cross(y))
|
||||
};
|
||||
Some(s)
|
||||
}
|
||||
|
||||
fn from_v3(v: mt::Vector3) -> Vec3 {
|
||||
Vec3::new(v.x, v.y, v.z)
|
||||
}
|
||||
|
||||
fn knots(k: &mt::KnotVec) -> Option<Vec<f64>> {
|
||||
let v: Vec<f64> = k.iter().copied().collect();
|
||||
(v.len() <= MAX_EVAL_KNOTS && v.iter().all(|x| x.is_finite())).then_some(v)
|
||||
}
|
||||
|
||||
/// Truck's `control_points[i][j]` (i along u) as the evaluator's u-fastest array.
|
||||
fn nurbs((pu, pv): (usize, usize), kv: &(mt::KnotVec, mt::KnotVec), rows: Vec<Vec<[f64; 4]>>) -> Option<EvalSurface> {
|
||||
let nu = rows.len();
|
||||
let nv = rows.first()?.len();
|
||||
if pu == 0 || pv == 0 || pu > MAX_SPLINE_DEGREE || pv > MAX_SPLINE_DEGREE || nu.checked_mul(nv)? > MAX_EVAL_POINTS {
|
||||
return None;
|
||||
}
|
||||
let mut ctrl = Vec::with_capacity(nu * nv);
|
||||
for j in 0..nv {
|
||||
for row in &rows {
|
||||
ctrl.push(<[[f64; 4]]>::get(row, j).copied()?);
|
||||
}
|
||||
}
|
||||
if ctrl.iter().flatten().any(|x| !x.is_finite()) {
|
||||
return None;
|
||||
}
|
||||
let knots = [knots(&kv.0)?, knots(&kv.1)?];
|
||||
if !(knots[0].len() == nu + pu + 1 && knots[1].len() == nv + pv + 1) {
|
||||
return None;
|
||||
}
|
||||
let closed = [closed_along(0, &ctrl, [nu, nv], [pu, pv], &knots), closed_along(1, &ctrl, [nu, nv], [pu, pv], &knots)];
|
||||
Some(EvalSurface::Nurbs { degree: [pu, pv], knots, count: [nu, nv], ctrl, closed })
|
||||
}
|
||||
|
||||
/// Does the surface join itself along direction `d` (a tube's way round)? Its knots there must
|
||||
/// be clamped, so the ends of the range are the first and last rows of control points, and
|
||||
/// those two rows must be the same points: the join is then the parameter line `d = end`, and
|
||||
/// the direction is periodic.
|
||||
fn closed_along(d: usize, ctrl: &[[f64; 4]], count: [usize; 2], degree: [usize; 2], knots: &[Vec<f64>; 2]) -> bool {
|
||||
let (n, p) = (count[d], degree[d]);
|
||||
if n < 2 {
|
||||
return false;
|
||||
}
|
||||
let Some(k) = knots.get(d) else { return false };
|
||||
let Some((&a, &b)) = k.first().zip(k.last()) else { return false };
|
||||
// Clamped: degree + 1 knots equal at each end, and a range to speak of.
|
||||
if !(b > a && k.iter().take(p + 1).all(|x| *x == a) && k.iter().rev().take(p + 1).all(|x| *x == b)) {
|
||||
return false;
|
||||
}
|
||||
// Control point `(i, j)` with `i` along u, as a point.
|
||||
let at =
|
||||
|i: usize, j: usize| ctrl.get(j * count[0] + i).and_then(|h| (h[3].abs() > 1e-300).then(|| Vec3::new(h[0] / h[3], h[1] / h[3], h[2] / h[3])));
|
||||
let along_d = |x: usize, t: usize| if d == 0 { at(x, t) } else { at(t, x) };
|
||||
let (mut extent, mut gap) = (0.0f64, 0.0f64);
|
||||
for t in 0..count[1 - d] {
|
||||
let (Some(first), Some(last)) = (along_d(0, t), along_d(n - 1, t)) else { return false };
|
||||
extent = extent.max(first.len()).max(last.len());
|
||||
gap = gap.max(first.dist(last));
|
||||
}
|
||||
gap <= 1e-9 * extent.max(1.0)
|
||||
}
|
||||
|
||||
/// A profile curve in the evaluator's form (a line extends for ever, as truck's does).
|
||||
pub(crate) fn eval_curve(c: &mt::Curve) -> Option<EvalCurve> {
|
||||
let (degree, kv, ctrl): (usize, Vec<f64>, Vec<[f64; 4]>) = match c {
|
||||
mt::Curve::Line(l) => (1, vec![0.0, 0.0, 1.0, 1.0], vec![[l.0.x, l.0.y, l.0.z, 1.0], [l.1.x, l.1.y, l.1.z, 1.0]]),
|
||||
mt::Curve::BSplineCurve(b) => (b.degree(), knots(b.knot_vec())?, b.control_points().iter().map(|p| [p.x, p.y, p.z, 1.0]).collect()),
|
||||
mt::Curve::NurbsCurve(n) => (n.degree(), knots(n.knot_vec())?, n.control_points().iter().map(|h| [h.x, h.y, h.z, h.w]).collect()),
|
||||
mt::Curve::IntersectionCurve(_) => return None,
|
||||
let revolved = |profile: EvalCurve, origin: Vec3, axis: Vec3| -> Result<EvalSurface> {
|
||||
let axis = axis.normalized().ok_or_else(|| invalid("an analytic surface's frame"))?;
|
||||
Ok(EvalSurface::Revolved { profile, origin, axis, transform: IDENTITY, swap: false })
|
||||
};
|
||||
let ok = (1..=MAX_SPLINE_DEGREE).contains(°ree)
|
||||
&& ctrl.len() <= MAX_EVAL_POINTS
|
||||
&& kv.len() == ctrl.len() + degree + 1
|
||||
&& ctrl.iter().flatten().all(|x| x.is_finite());
|
||||
ok.then_some(EvalCurve { degree, knots: kv, ctrl })
|
||||
// The generator of a cylinder or cone over the face's heights, run from the top down so the
|
||||
// revolution's (generator, angle) normal is the record's natural normal (outward).
|
||||
let generator = |at: &dyn Fn(f64) -> Vec3| -> Result<EvalCurve> {
|
||||
let [h0, h1] = domain[1];
|
||||
if !(h0.is_finite() && h1.is_finite() && h1 > h0) {
|
||||
return Err(invalid("a cylinder or cone face with no height"));
|
||||
}
|
||||
let (a, b) = (at(h1), at(h0));
|
||||
// Parameter −v: the generator at t is the surface at height −t.
|
||||
Ok(EvalCurve { degree: 1, knots: vec![-h1, -h1, -h0, -h0], ctrl: vec![[a.x, a.y, a.z, 1.0], [b.x, b.y, b.z, 1.0]] })
|
||||
};
|
||||
match surface {
|
||||
brep::Surface::Cylinder { frame, radius } => {
|
||||
let (o, x, z) = axes(frame);
|
||||
revolved(generator(&|h| o + x * *radius + z * h)?, o, z)
|
||||
}
|
||||
brep::Surface::Cone { frame, radius_at_origin, half_angle } => {
|
||||
let (o, x, z) = axes(frame);
|
||||
let slope = half_angle.tan();
|
||||
revolved(generator(&|h| o + x * (radius_at_origin + h * slope) + z * h)?, o, z)
|
||||
}
|
||||
brep::Surface::Sphere { frame, radius } => {
|
||||
let (o, x, z) = axes(frame);
|
||||
// The meridian from the north pole to the south pole.
|
||||
revolved(arc(o, x, z, *radius, [-FRAC_PI_2, FRAC_PI_2]).reversed(), o, z)
|
||||
}
|
||||
brep::Surface::Torus { frame, major_radius, minor_radius } => {
|
||||
let (o, x, z) = axes(frame);
|
||||
revolved(arc(o + x * *major_radius, x, z, *minor_radius, [0.0, TAU]).reversed(), o, z)
|
||||
}
|
||||
_ => match SurfaceEval::new(surface) {
|
||||
Ok(SurfaceEval::Eval(e, _)) => Ok(e),
|
||||
Ok(_) => Err(invalid("an analytic surface out of place")),
|
||||
Err(e) => Err(KernelError::Failed(format!("not supported yet: tessellation: {e}"))),
|
||||
},
|
||||
}
|
||||
}
|
||||
|
||||
/// The circle `centre + r (cos t x + sin t y)` for t in `angle`, as an exact rational quadratic
|
||||
/// spline over [0, 1], one span per quarter turn or less.
|
||||
fn arc(centre: Vec3, x: Vec3, y: Vec3, r: f64, angle: [f64; 2]) -> EvalCurve {
|
||||
let span = angle[1] - angle[0];
|
||||
let count = (span / std::f64::consts::FRAC_PI_2).ceil().clamp(1.0, 4.0) as usize;
|
||||
let step = span / count as f64;
|
||||
let weight = (step / 2.0).cos();
|
||||
// A point at angle `t`, its radius divided by `scale` (the middle control point of a span
|
||||
// lies 1 / cos(step / 2) out), as homogeneous coordinates of weight `w`.
|
||||
let at = |t: f64, scale: f64, w: f64| {
|
||||
let p = centre + x * (r * t.cos() / scale) + y * (r * t.sin() / scale);
|
||||
[p.x * w, p.y * w, p.z * w, w]
|
||||
};
|
||||
let mut ctrl = vec![at(angle[0], 1.0, 1.0)];
|
||||
let mut knots = vec![0.0; 3];
|
||||
for i in 0..count {
|
||||
let start = angle[0] + i as f64 * step;
|
||||
ctrl.push(at(start + step / 2.0, weight, weight));
|
||||
ctrl.push(at(start + step, 1.0, 1.0));
|
||||
knots.extend(std::iter::repeat_n((i + 1) as f64 / count as f64, if i + 1 == count { 3 } else { 2 }));
|
||||
}
|
||||
EvalCurve { degree: 2, knots, ctrl }
|
||||
}
|
||||
|
||||
/// A curve on [0, 1] run the other way: t becomes 1 − t.
|
||||
trait Reversed {
|
||||
fn reversed(self) -> Self;
|
||||
}
|
||||
|
||||
impl Reversed for EvalCurve {
|
||||
fn reversed(mut self) -> Self {
|
||||
self.ctrl.reverse();
|
||||
self.knots = self.knots.iter().rev().map(|k| 1.0 - k).collect();
|
||||
self
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
|
||||
Reference in New Issue
Block a user