kernel: a face's boundary takes its (u, v) from the record's exact pcurves

Each edge node keeps its curve parameter; each face reads its nodes' parameters from its own
pcurve of the edge at that parameter, so nothing on the boundary is projected and a seam's two
uses (each with its pcurve) are its two sides: the cases the projection could not tell apart
("boundary winds around the surface", sides that cross where projection snapped them, poles).

The record's (u, v) go to the evaluator's in closed form (tess::surface::Params): the same for
planes, splines and revolutions; (θ, v) to (−v, θ) on a cylinder's or cone's generator; on a
sphere or torus (θ, φ) to (1 − s(φ), θ), s the rational meridian's parameter, from
tan(φ'/2) running linearly over each span (φ' from the span's middle).

A pcurve that does not cover its edge's interval, or places a node farther from the edge than
ten times the chord tolerance plus the edge's own tolerance, is an error naming the face and the
edge: the record is wrong there (for an analytic face, a trim left in SMLib's internal NURBS
parameters rather than the record's would show here).

(not built: base in progress) tess/ has no errors under cargo check.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_017itGh5fpMEousnZ91cYcVt
This commit is contained in:
jyviko
2026-10-22 00:32:36 -04:00
parent f73b3be1a4
commit d61efb306c
3 changed files with 177 additions and 115 deletions
+15 -10
View File
@@ -16,6 +16,9 @@ const MAX_DEPTH: u32 = 24;
/// An edge's nodes, from its start vertex (interval min) to its end vertex (interval max).
pub(crate) struct Nodes {
pub points: Vec<Vec3>,
/// The curve's parameter at each node (the edge's interval, ends included): where each
/// face's pcurve of the edge gives the node's (u, v).
pub params: Vec<f64>,
}
/// Nodes of edge `ei` of `b` from its start to its end: every chord within the tolerance of
@@ -32,6 +35,10 @@ pub(crate) fn sample_edge(b: &brep::Brep, ei: usize, rules: &Rules, h: f64) -> R
};
let (front, back) = (vertex(0)?, vertex(1)?);
let tol = rules.tol;
let [t0, t1] = edge.interval;
if !(t0.is_finite() && t1.is_finite() && t1 > t0) {
return Err(fail(format!("interval [{t0}, {t1}]")));
}
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
@@ -39,14 +46,12 @@ pub(crate) fn sample_edge(b: &brep::Brep, ei: usize, rules: &Rules, h: f64) -> R
// forwards takes the next edge's first node (the vertex) — a crack.
let n = if h > 0.0 { (front.dist(back) / h).ceil().clamp(1.0, (MAX_EDGE_NODES - 1) as f64) as usize } else { 1 };
let mut pts: Vec<Vec3> = (0..=n).map(|i| front + (back - front) * (i as f64 / n as f64)).collect();
if let Some(p) = pts.last_mut() {
let mut params: Vec<f64> = (0..=n).map(|i| t0 + (t1 - t0) * (i as f64 / n as f64)).collect();
if let (Some(p), Some(t)) = (pts.last_mut(), params.last_mut()) {
*p = back;
*t = t1;
}
return Ok(Nodes { points: pts });
}
let [t0, t1] = edge.interval;
if !(t0.is_finite() && t1.is_finite() && t1 > t0) {
return Err(fail(format!("interval [{t0}, {t1}]")));
return Ok(Nodes { points: pts, params });
}
let curve = CurveEval::new(&edge.curve).map_err(|e| fail(e.to_string()))?;
let at =
@@ -67,7 +72,7 @@ pub(crate) fn sample_edge(b: &brep::Brep, ei: usize, rules: &Rules, h: f64) -> R
mid: Option<Vec3>,
depth: u32,
}
let mut out: Vec<Vec3> = vec![front];
let mut out: Vec<(f64, Vec3)> = vec![(t0, front)];
let mut stack: Vec<Seg> = Vec::with_capacity(64);
for (i, w) in ends.windows(2).enumerate().rev() {
let (Some(&pa), Some(&pb)) = (w.first(), w.get(1)) else { continue };
@@ -94,13 +99,13 @@ pub(crate) fn sample_edge(b: &brep::Brep, ei: usize, rules: &Rules, h: f64) -> R
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);
out.push((b, pb));
}
}
if let Some(p) = out.last_mut() {
*p = back;
*p = (t1, back);
}
Ok(Nodes { points: out })
Ok(Nodes { points: out.iter().map(|p| p.1).collect(), params: out.iter().map(|p| p.0).collect() })
}
/// Angle (radians) between two directions; zero when either is empty.
+87 -86
View File
@@ -1,6 +1,8 @@
//! 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.
//! exactly, handed over as a [`TessPlan`] for the CPU or GPU evaluator. A boundary node's
//! parameters on each face come from that face's exact pcurve of the edge, at the node's curve
//! parameter: nothing is projected, and a seam's two sides are its two pcurves.
//!
//! Per face:
//! 1. a cheaper layout where the face suits one: planar faces whose boundary runs along their
@@ -12,15 +14,15 @@
//! 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.
//! seam or pole), or whose pcurves do not lie on its edges, is an error naming the face (and
//! the edge) and why. Its cause is upstream, in the record's pcurves or seams, and is fixed
//! there; there is no other mesher to fall back to.
pub(crate) mod cdt;
pub(crate) mod collar;
pub(crate) mod coons;
pub(crate) mod cutgrid;
pub(crate) mod edges;
pub(crate) mod fallback;
pub(crate) mod grid;
pub(crate) mod motorcycle;
pub(crate) mod pieces;
@@ -44,8 +46,6 @@ pub struct PlanStats {
pub planar: usize,
pub structured: usize,
pub delaunay: usize,
/// Faces no method could plan, filled flat from their boundary nodes.
pub filled: usize,
}
/// Curved faces from which their sizes are worked out on the worker pool.
@@ -64,15 +64,6 @@ const PARALLEL_EDGES: usize = 256;
const MAX_LINES: usize = 8192;
/// Most points Ruppert refinement may add to one face, beyond its seeds.
const MAX_REFINE: usize = 20_000;
/// Two boundary nodes this close together on the surface (mm) are one node. Far below any
/// spacing a mesh can carry, and far above the parameters a projection tells apart. See
/// [`drop_collapsed`].
const SAME_NODE: f64 = 1e-6;
/// How far a face's boundary may lie off its surface and still be planned on it, as a fraction of
/// the body's extent, at every display quality (ten chord tolerances at Normal).
const ON_SURFACE: f64 = 1e-2;
/// How much finer a refused face's edges are divided, retry by retry.
const REFINE_STEPS: [f64; 2] = [3.0, 10.0];
fn unsupported(what: impl std::fmt::Display) -> KernelError {
KernelError::Failed(format!("not supported yet: tessellation: {what}"))
@@ -141,11 +132,12 @@ pub(crate) fn plan(b: &brep::Brep, tol: f64, turn: Option<size::Turn>) -> Result
// 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<EvalSurface> = faces
let evals: Vec<(EvalSurface, surface::Params)> = 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 (evals, params): (Vec<EvalSurface>, Vec<surface::Params>) = evals.into_iter().unzip();
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)))
@@ -163,46 +155,12 @@ pub(crate) fn plan(b: &brep::Brep, tol: f64, turn: Option<size::Turn>) -> 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 {
let (Some(face), Some(eval), Some(map), Some(size)) = (faces.get(fi), evals.get(fi), params.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();
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();
}
ep.pop();
pts.extend(ep);
}
if pts.len() >= 2 {
loops3.push(pts);
}
}
Ok((loops3, known))
};
// A plan in the surface's orientation turned the face's way.
let oriented_face = |face: &mt::Face, fp: FacePlan| FacePlan {
flip: !face.orientation(),
triangles: if face.orientation() { fp.triangles } else { fp.triangles.into_iter().map(|[a, b, c]| [a, c, b]).collect() },
..fp
};
let surface_of = |fi: usize| match sized.get(fi) {
Some((Some(eval), Some(size))) => Ok((eval.clone(), *size)),
_ => Err(unsupported(format!("face {fi}: surface not evaluable"))),
};
// One face, from the edge nodes as they are and with the rules for its level of detail.
let plan_one = |fi: usize, face: &mt::Face, nodes: &HashMap<mt::EdgeID, edges::Nodes>, rules: &size::Rules| {
let (eval, size) = surface_of(fi)?;
let (loops3, known) = boundary(face, &eval, nodes)?;
let loops = boundary(b, fi, face, eval, map, &nodes, &rules)?;
let mut stats = PlanStats::default();
let fp = plan_face(fi, kind(&face.surface), eval, &loops3, &rules, *size, &mut stats)?;
let fp = plan_face(fi, kind(&face.surface), eval, &loops, &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 };
@@ -218,13 +176,83 @@ pub(crate) fn plan(b: &brep::Brep, tol: f64, turn: Option<size::Turn>) -> Result
stats.planar += s.planar;
stats.structured += s.structured;
stats.delaunay += s.delaunay;
stats.filled += s.filled;
}
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))
}
/// Face `fi`'s boundary, loop by loop: the shared edge nodes (each edge's last node is the next
/// edge's first) and their parameters on the face's surface, from each edge use's exact pcurve
/// at the node's curve parameter. A seam edge's two uses on one face have their own pcurves, so
/// its two sides are told apart. Each loop is unwrapped along its periodic directions (a pcurve
/// may start a period along from where the previous one ended). A vertex loop (a pole) bounds no
/// area and has no nodes of its own. A pcurve that does not cover its edge's interval, or puts a
/// node farther off the edge than the tolerances allow, is an error naming the face and edge:
/// the record is wrong there.
fn boundary(
b: &brep::Brep,
fi: usize,
face: &brep::Face,
s: &EvalSurface,
map: &surface::Params,
nodes: &[edges::Nodes],
rules: &size::Rules,
) -> Result<Vec<seam::Loop>> {
let dom = uv::domain(s);
let fail = |e: u32, why: String| unsupported(format!("face {fi} ({}): edge {e}: {why}", kind(&face.surface)));
let mut loops = Vec::with_capacity(face.loops.len());
for l in &face.loops {
let Loop::Edges { edges } = l else { continue };
let (mut pts, mut uvs): (Vec<Vec3>, Vec<[f64; 2]>) = (Vec::new(), Vec::new());
for oe in edges {
let n = nodes.get(oe.edge as usize).ok_or_else(|| fail(oe.edge, "no nodes".into()))?;
let edge = b.edges.get(oe.edge as usize).ok_or_else(|| fail(oe.edge, "no such edge".into()))?;
let [t0, t1] = edge.interval;
let [p0, p1] = oe.pcurve.interval().ok_or_else(|| fail(oe.edge, "a malformed pcurve".into()))?;
let slack = 1e-9 * (t1 - t0).abs().max(t0.abs()).max(t1.abs()).max(1.0);
if (p0 - t0).abs() > slack || (p1 - t1).abs() > slack {
return Err(fail(oe.edge, format!("its pcurve spans [{p0}, {p1}], the edge [{t0}, {t1}]")));
}
// The edge's own tolerance: how far its curve may stand off the face.
let accept = (rules.tol * 10.0).max(1e-6) + edge.tolerance.max(0.0);
let mut ep: Vec<Vec3> = n.points.clone();
let mut eq: Vec<[f64; 2]> = Vec::with_capacity(ep.len());
for (k, (t, p)) in n.params.iter().zip(&n.points).enumerate() {
let q = oe.pcurve.eval(t.clamp(p0, p1)).ok_or_else(|| fail(oe.edge, format!("its pcurve does not evaluate at t = {t}")))?;
let q = map.to_eval(q);
let on =
s.derivs(q[0], q[1]).map(|d| d.0).ok_or_else(|| fail(oe.edge, format!("node {k}: the surface does not evaluate at {q:?}")))?;
if !(on.dist(*p) <= accept) {
return Err(fail(oe.edge, format!("node {k}: its pcurve is {} mm off the edge at t = {t}", on.dist(*p))));
}
eq.push(q);
}
if !oe.forward {
ep.reverse();
eq.reverse();
}
ep.pop();
eq.pop();
pts.extend(ep);
uvs.extend(eq);
}
let mut prev: Option<[f64; 2]> = None;
for q in &mut uvs {
if let Some(pr) = prev {
for d in 0..2 {
q[d] = uv::unwrap(q[d], pr[d], dom.period[d]);
}
}
prev = Some(*q);
}
if pts.len() >= 2 {
loops.push(seam::Loop { uv: uvs, pts });
}
}
Ok(loops)
}
/// 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
@@ -413,7 +441,7 @@ fn plan_face(
fi: usize,
kind: &str,
s: &EvalSurface,
loops3: &[Vec<Vec3>],
loops: &[seam::Loop],
rules: &size::Rules,
size: (f64, [f64; 2]),
stats: &mut PlanStats,
@@ -421,35 +449,8 @@ fn plan_face(
let tol = rules.tol;
let (h, speed) = size;
let dom = uv::domain(s);
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() {
// 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]);
}
}
luv.push(q);
prev = Some(q);
}
let mut l = seam::Loop { uv: luv, pts: pts.clone() };
drop_collapsed(&mut l, speed);
for (d, per) in dom.period.iter().enumerate() {
if let Some(per) = per {
seam::close_through_pole(s, &mut l, d, *per);
}
}
projected.push(l);
}
align_periods(&mut projected, &dom);
let mut projected: Vec<seam::Loop> = loops.to_vec();
// A loop that winds around a periodic direction does not close in (u, v): the face is cut
// along a parameter line into one period (a seamless cylinder face, a helical trim).
// A loop that runs up a seam from a pole and straight back down it (a half sphere bounded
@@ -511,7 +512,7 @@ fn plan_face(
}
ranges.push(start..uvs.len());
}
if uvs.len() < MIN_BOUNDARY_NODES {
if uvs.len() < 3 {
return Err(fail("degenerate boundary"));
}
let mut bbox = [(f64::MAX, f64::MIN); 2];
@@ -529,7 +530,7 @@ fn plan_face(
let seeds = cdt::lattice(bbox, step, &polys, speed);
let eval = |u: f64, v: f64| s.derivs(u, v).map(|(p, su, sv)| (p, su.cross(sv).normalized()));
let sag = cdt::Sag { eval: &eval, tol, h, turn: rules.turn, scale: speed };
let fill = cdt::Fill { seeds: &seeds, refine: Some((sag, MAX_REFINE)), shape: !matches!(s, EvalSurface::Plane { .. }), lenient: false };
let fill = cdt::Fill { seeds: &seeds, refine: Some((sag, MAX_REFINE)), shape: !matches!(s, EvalSurface::Plane { .. }) };
cdt::triangulate(&uvs, &ranges, fill).ok_or_else(|| fail("boundary crosses itself"))
};
// Planar and curved faces alike: a structured grid where the boundary runs along the
@@ -766,7 +767,7 @@ fn plan_pieces(f: &Region) -> Option<(Vec<[f64; 2]>, Vec<[u32; 3]>)> {
let seeds = cdt::lattice(bbox, [f.h / f.speed[0], f.h / f.speed[1]], &polys, f.speed);
let eval = |u: f64, v: f64| s.derivs(u, v).map(|(p, su, sv)| (p, su.cross(sv).normalized()));
let sag = cdt::Sag { eval: &eval, tol: f.rules.tol, h: f.h, turn: f.rules.turn, scale: f.speed };
cdt::triangulate(&luv, &ranges, cdt::Fill { seeds: &seeds, refine: Some((sag, MAX_REFINE)), shape: true, lenient: false })
cdt::triangulate(&luv, &ranges, cdt::Fill { seeds: &seeds, refine: Some((sag, MAX_REFINE)), shape: true })
})?;
let (lextra, ltris) = local;
let base = u32::try_from(n + extra.len()).ok()?;
+75 -19
View File
@@ -1,10 +1,11 @@
//! 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.
//! planner and the CPU and GPU evaluators share, and the map from the record's (u, v) to the
//! evaluator's, so a face's boundary takes its parameters from its edges' exact pcurves.
//! 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 by the
//! record's angle u; the generator's parameter is −v on a line and, on an arc, the rational
//! parameter of the record's angle v (closed form, [`Arc::param`]).
use solidcraft_geom::Vec3;
use solidcraft_geom::brep::eval::SurfaceEval;
@@ -20,10 +21,59 @@ fn v3(p: [f64; 3]) -> Vec3 {
/// 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> {
/// How a face's record parameters become its evaluator's.
#[derive(Clone, Copy, Debug, PartialEq)]
pub(crate) enum Params {
/// The same parameters (planes, splines, revolutions).
Same,
/// A cylinder or cone: (angle, height) to (−height, angle).
Height,
/// A sphere or torus: (angle, meridian angle) to (the reversed arc's parameter, angle).
Meridian(Arc),
}
impl Params {
/// The evaluator's parameters of the record's `uv`.
pub(crate) fn to_eval(&self, uv: [f64; 2]) -> [f64; 2] {
let [u, v] = uv;
match self {
Params::Same => uv,
Params::Height => [-v, u],
Params::Meridian(a) => [1.0 - a.param(v), u],
}
}
}
/// The spans of a circular arc made by [`arc`]: from angle `start`, `count` spans of `step`
/// each, over the parameter [0, 1]; `periodic` for a whole turn.
#[derive(Clone, Copy, Debug, PartialEq)]
pub(crate) struct Arc {
start: f64,
step: f64,
count: usize,
periodic: bool,
}
impl Arc {
/// The arc's parameter at `angle`. On a rational quadratic span of half angle α with middle
/// weight cos α, tan(θ / 2) runs linearly from −tan(α / 2) to tan(α / 2) with the span's
/// own parameter (θ from the span's middle). A whole turn continues by whole parameters.
pub(crate) fn param(&self, angle: f64) -> f64 {
let total = self.step * self.count as f64;
let turns = if self.periodic { ((angle - self.start) / total).floor() } else { 0.0 };
let a = angle - turns * total;
let i = (((a - self.start) / self.step).floor().max(0.0) as usize).min(self.count.saturating_sub(1));
let half = 0.5 * self.step;
let local = a - (self.start + (i as f64 + 0.5) * self.step);
let s = 0.5 * ((0.5 * local).tan() / (0.5 * half).tan() + 1.0);
(i as f64 + s) / self.count as f64 + turns
}
}
/// `surface` in the evaluator's form, and how the record's parameters map onto it. `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.
pub(crate) fn eval_surface(surface: &brep::Surface, domain: [[f64; 2]; 2]) -> Result<(EvalSurface, Params)> {
use std::f64::consts::{FRAC_PI_2, TAU};
let invalid = |what: &str| KernelError::Invalid(format!("tessellation: {what}"));
let axes = |f: &Frame| {
@@ -45,27 +95,32 @@ pub(crate) fn eval_surface(surface: &brep::Surface, domain: [[f64; 2]; 2]) -> Re
// 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]] })
};
// A meridian from the record's least angle to its greatest, reversed (so the normal is the
// record's natural one, outward).
let meridian = |centre: Vec3, x: Vec3, z: Vec3, r: f64, angle: [f64; 2], o: Vec3| -> Result<(EvalSurface, Params)> {
let (curve, spans) = arc(centre, x, z, r, angle);
Ok((revolved(curve.reversed(), o, z)?, Params::Meridian(spans)))
};
match surface {
brep::Surface::Cylinder { frame, radius } => {
let (o, x, z) = axes(frame);
revolved(generator(&|h| o + x * *radius + z * h)?, o, z)
Ok((revolved(generator(&|h| o + x * *radius + z * h)?, o, z)?, Params::Height))
}
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)
Ok((revolved(generator(&|h| o + x * (radius_at_origin + h * slope) + z * h)?, o, z)?, Params::Height))
}
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)
meridian(o, x, z, *radius, [-FRAC_PI_2, FRAC_PI_2], o)
}
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)
meridian(o + x * *major_radius, x, z, *minor_radius, [0.0, TAU], o)
}
_ => match SurfaceEval::new(surface) {
Ok(SurfaceEval::Eval(e, _)) => Ok(e),
Ok(SurfaceEval::Eval(e, _)) => Ok((e, Params::Same)),
Ok(_) => Err(invalid("an analytic surface out of place")),
Err(e) => Err(KernelError::Failed(format!("not supported yet: tessellation: {e}"))),
},
@@ -73,8 +128,8 @@ pub(crate) fn eval_surface(surface: &brep::Surface, domain: [[f64; 2]; 2]) -> Re
}
/// 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 {
/// spline over [0, 1], one span per quarter turn or less, and its spans.
fn arc(centre: Vec3, x: Vec3, y: Vec3, r: f64, angle: [f64; 2]) -> (EvalCurve, Arc) {
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;
@@ -93,7 +148,8 @@ fn arc(centre: Vec3, x: Vec3, y: Vec3, r: f64, angle: [f64; 2]) -> EvalCurve {
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 }
let periodic = (span - std::f64::consts::TAU).abs() < 1e-12;
(EvalCurve { degree: 2, knots, ctrl }, Arc { start: angle[0], step, count, periodic })
}
/// A curve on [0, 1] run the other way: t becomes 1 − t.