Skip to main content

brep_kernel/blending/blend/chain/
closed.rs

1use crate::{KernelRefusal, KernelStage, OrRefuse};
2use super::*;
3
4/// Piece of a support rim assigned to one mate face.
5pub(super) struct RimPiece {
6    pub(super) face_id: u64,
7    pub(super) loop_index: usize,
8    /// Global-parameter window of this piece.
9    pub(super) window: [f64; 2],
10    pub(super) edge_id: u64,
11}
12
13/// Blend a chain of conjugated edges with one rolling-ball blend face
14/// (Golovanov §6.9.5: all conjugated edges processed together).  A CLOSED
15/// chain (stadium rim, T-pipe saddle) welds into a periodic blend; an OPEN
16/// chain (line->arc->line capped by end faces) rolls a clamped blend that
17/// terminates in a transverse edge on each end face.
18pub fn blend_smooth_chain(
19    solid: &BrepSolid,
20    seed_edge_id: u64,
21    radius: f64,
22    chamfer: bool,
23    name: Option<&str>,
24) -> Result<BrepSolid, KernelRefusal> {
25    // One blend operation: the reports of an earlier one, never consumed, are
26    // dropped as it starts (nested calls keep this operation's).
27    let _operation = crate::blend::BlendOperation::enter();
28    if !(radius > 0.0) || !radius.is_finite() {
29        return Err(KernelRefusal::input(KernelStage::Collect, "radius", "blend: radius must be positive"));
30    }
31    let chain = collect_smooth_chain(solid, seed_edge_id)?;
32    if chain.segments.len() < 2 {
33        return Err(KernelRefusal::internal(KernelStage::Refine, "chain_single_segment", "blend: chain collapsed to a single segment"));
34    }
35    if chain.closed {
36        blend_closed_smooth_chain(solid, &chain.segments, radius, chamfer, name)
37    } else {
38        blend_open_smooth_chain(solid, &chain, radius, chamfer, name)
39    }
40}
41
42/// `blend_smooth_chain` restricted to CLOSED chains: an edge that is one arc
43/// of a seam-split rim (a cyl×cyl saddle) is blended as the whole rim, but an
44/// OPEN tangent chain is not walked — filleting one edge of it must not
45/// fillet the unselected edges it runs into.  Tangent propagation is a
46/// selection policy for the feature layer; the kernel caps an unselected
47/// continuation instead (`network.rs`).
48pub fn blend_smooth_chain_if_closed(
49    solid: &BrepSolid,
50    seed_edge_id: u64,
51    radius: f64,
52    chamfer: bool,
53    name: Option<&str>,
54) -> Result<BrepSolid, KernelRefusal> {
55    // One blend operation: the reports of an earlier one, never consumed, are
56    // dropped as it starts (nested calls keep this operation's).
57    let _operation = crate::blend::BlendOperation::enter();
58    if !(radius > 0.0) || !radius.is_finite() {
59        return Err(KernelRefusal::input(KernelStage::Collect, "radius", "blend: radius must be positive"));
60    }
61    let chain = collect_smooth_chain(solid, seed_edge_id)?;
62    if chain.segments.len() < 2 {
63        return Err(KernelRefusal::internal(KernelStage::Refine, "chain_single_segment", "blend: chain collapsed to a single segment"));
64    }
65    if !chain.closed {
66        return Err(
67            KernelRefusal::unsupported(KernelStage::Classify, "open_tangent_chain", "blend: the edge is one arc of an OPEN tangent chain; the unselected \
68             continuation is capped, not blended"),
69        );
70    }
71    blend_closed_smooth_chain(solid, &chain.segments, radius, chamfer, name)
72}
73
74// BREP lane L: begin lcount236_chain_ladder_memo_20261016
75/// A closed chain's ladder that tops out answers the same refusal for the same
76/// inputs, and a fillet selection asks it again and again: the group's cutter
77/// composition, then every drop-one subset, each blending the same collar on a
78/// solid the earlier edges did not touch there (up to four calls a case). Up
79/// to its verdict the ladder reads only the chain — every segment's edge, both
80/// mate faces (surface, sense, loops) and coedges, the signed radii — the
81/// radius, the chamfer flag, the bar and three switches, never the rest of the
82/// solid; the key is exactly those, bit for bit. Only the `chain_ladder_top`
83/// refusal is kept, and a hit returns that refusal unchanged.
84/// `BREP_CHAIN_LADDER_MEMO=0` turns the memo off; `=verify` recomputes every
85/// hit and panics unless the recomputed answer is the cached refusal, byte for
86/// byte.
87fn blend_closed_smooth_chain(
88    solid: &BrepSolid,
89    segments: &[ChainSegment<'_>],
90    radius: f64,
91    chamfer: bool,
92    name: Option<&str>,
93) -> Result<BrepSolid, KernelRefusal> {
94    let mode = std::env::var("BREP_CHAIN_LADDER_MEMO");
95    let (memo_on, verify) = match mode.as_deref() {
96        Ok("0") => (false, false),
97        Ok("verify") => (true, true),
98        _ => (true, false),
99    };
100    if !memo_on {
101        return blend_closed_smooth_chain_unmemoized(solid, segments, radius, chamfer, name);
102    }
103    let bar = crate::KernelTolerances::for_solid(solid, 1e-7).intersection_fit;
104    let key = ladder_memo_key(segments, radius, chamfer, bar);
105    let cached = LADDER_MEMO.with(|memo| memo.borrow().iter().find(|(stored, _)| *stored == key).map(|(_, refusal)| refusal.clone()));
106    if let (Some(refusal), false) = (&cached, verify) {
107        crate::blend::record_count("count.chain_ladder_memo.hit", 1);
108        return Err(refusal.clone());
109    }
110    let result = blend_closed_smooth_chain_unmemoized(solid, segments, radius, chamfer, name);
111    if let Some(refusal) = cached {
112        let same = matches!(&result, Err(again) if *again == refusal && format!("{again:?}") == format!("{refusal:?}"));
113        assert!(same, "chain ladder memo mismatch: cached {refusal:?}, recomputed {:?}", result.as_ref().err());
114        crate::blend::record_count("count.chain_ladder_memo.verified", 1);
115        return result;
116    }
117    if let Err(refusal) = &result {
118        if matches!(&refusal.class, crate::RefusalClass::NonConvergence { what } if what == "chain_ladder_top") {
119            LADDER_MEMO.with(|memo| {
120                let mut memo = memo.borrow_mut();
121                if memo.len() >= LADDER_MEMO_CAPACITY {
122                    memo.remove(0);
123                }
124                memo.push((key, refusal.clone()));
125            });
126        }
127    }
128    result
129}
130
131/// How many topped-out ladders the memo keeps (oldest dropped first).
132const LADDER_MEMO_CAPACITY: usize = 32;
133
134thread_local! {
135    /// Topped-out closed-chain ladders: the bit-exact key and the refusal.
136    static LADDER_MEMO: std::cell::RefCell<Vec<(Vec<u64>, KernelRefusal)>> = const { std::cell::RefCell::new(Vec::new()) };
137}
138
139/// Every input the closed-chain ladder reads up to its verdict, as raw bits.
140fn ladder_memo_key(segments: &[ChainSegment<'_>], radius: f64, chamfer: bool, bar: f64) -> Vec<u64> {
141    fn curve(key: &mut Vec<u64>, curve: &NurbsCurve) {
142        key.extend([curve.degree as u64, curve.knots.len() as u64, curve.control_points.len() as u64]);
143        key.extend(curve.knots.iter().map(|knot| knot.to_bits()));
144        for point in &curve.control_points {
145            key.extend([point.x.to_bits(), point.y.to_bits(), point.z.to_bits(), point.w.to_bits()]);
146        }
147    }
148    fn face(key: &mut Vec<u64>, face: &FaceRecord) {
149        let surface = &face.surface;
150        key.extend([face.id, face.same_sense as u64, surface.degree_u as u64, surface.degree_v as u64]);
151        key.extend([surface.knots_u.len() as u64, surface.knots_v.len() as u64, surface.control_points.len() as u64]);
152        key.extend(surface.knots_u.iter().chain(&surface.knots_v).map(|knot| knot.to_bits()));
153        for row in &surface.control_points {
154            key.push(row.len() as u64);
155            for point in row {
156                key.extend([point.x.to_bits(), point.y.to_bits(), point.z.to_bits(), point.w.to_bits()]);
157            }
158        }
159        key.push(face.loops.len() as u64);
160        for loop_record in &face.loops {
161            key.extend([loop_record.id, loop_record.coedges.len() as u64]);
162            for coedge in &loop_record.coedges {
163                key.extend([coedge.id, coedge.edge_id, coedge.forward as u64]);
164                curve(key, &coedge.pcurve);
165            }
166        }
167    }
168    let mut key = vec![
169        radius.to_bits(),
170        chamfer as u64,
171        bar.to_bits(),
172        crate::blend::station_refinement_on() as u64,
173        crate::blend::chain_arc_stations_on() as u64,
174        wall_verdict_on() as u64,
175        local_continue_on() as u64,
176        segments.len() as u64,
177    ];
178    for segment in segments {
179        let edge = segment.edge;
180        key.extend([edge.id, edge.start_vertex_id, edge.end_vertex_id, edge.t0.to_bits(), edge.t1.to_bits(), segment.forward as u64]);
181        curve(&mut key, &edge.curve);
182        for mate in [&segment.first, &segment.second] {
183            key.extend([mate.coedge.id, mate.coedge.edge_id, mate.coedge.forward as u64, mate.loop_index as u64, mate.rho.to_bits()]);
184            curve(&mut key, &mate.coedge.pcurve);
185            face(&mut key, mate.face);
186        }
187    }
188    key
189}
190// BREP lane L: end lcount236_chain_ladder_memo_20261016
191
192/// Blend a CLOSED chain of conjugated edges with one rolling-ball blend
193/// face (Golovanov §6.9.5: all conjugated edges processed together).
194///
195/// A wall whose ball-centre curve turns tighter than the ball FOLDS, and the
196/// envelope a plain march sweeps is then not the blend (`blend/fold.rs`).  It
197/// is not refused here any more where the fold is a lens the wall closes on
198/// itself: the chain is re-marched up a measured ladder of budgets and the lens
199/// is carved out along the crease (`blend/carve.rs`).  The re-march is the
200/// whole of the extra cost and only a folding wall pays it — the first march is
201/// the ordinary one, and its refusal is what asks for the second.
202fn blend_closed_smooth_chain_unmemoized(
203    solid: &BrepSolid,
204    segments: &[ChainSegment<'_>],
205    radius: f64,
206    chamfer: bool,
207    name: Option<&str>,
208) -> Result<BrepSolid, KernelRefusal> {
209    let bar = crate::KernelTolerances::for_solid(solid, 1e-7).intersection_fit;
210    // The plain collar climbs the same measured ladder the carve does, from the
211    // same first rung: a budget picked once is a resolution request in disguise.
212    let mut budget = ChainBudget::for_chain(segments, CHAIN_PER_SEGMENT)?;
213    // Every too-coarse rung's reading, for the refusal at the top.
214    let mut trend: Vec<(usize, f64)> = Vec::new();
215    let first = loop {
216        let per_segment = budget.rung;
217        let rung = build_closed_smooth_chain(
218            solid,
219            segments,
220            radius,
221            chamfer,
222            name,
223            &budget,
224            FoldPolicy::Refuse,
225            bar,
226        );
227        match rung {
228            Ok(Rung::TooCoarse { deviation }) => {
229                crate::blend::carve::carve_trace(format_args!(
230                    "blend chain: {} leaves the rails \
231                     {deviation:.6e} from the rolling ball, against {bar:.3e}",
232                    budget.describe()
233                ));
234                trend.push((per_segment, deviation));
235                if budget.largest() * 2 > CHAIN_CARVE_MAX_PER_SEGMENT {
236                    let (best_rung, best, listing) = ladder_trend(&trend);
237                    return Err(KernelRefusal::non_convergence(KernelStage::Refine, "chain_ladder_top", format!(
238                        "blend: the chain's ladder reached its top, {per_segment} stations per segment, \
239                         and no rung passed — its best, {best_rung} stations per segment, left this \
240                         collar's fitted rails {best:.6e} from the rolling ball's own contacts against \
241                         a bar of {bar:.3e} (rung by rung: {listing}), so there is no wall accurate \
242                         enough to build"
243                    )));
244                }
245                budget = budget.doubled();
246            }
247            // The continuation Newton of one station did not converge from its
248            // neighbour at this spacing. Only that typed station failure climbs:
249            // the next rung seeds every station from a neighbour half as far
250            // away. At the top the march's own refusal is returned verbatim.
251            Ok(Rung::StationFailed(failure)) => {
252                crate::blend::carve::carve_trace(format_args!(
253                    "blend chain: {per_segment} stations per segment: segment {} position {} \
254                     (t {:.9}) continued from position {} (t {:.9}, seed {:?}) did not converge: {}",
255                    failure.segment, failure.position, failure.t, failure.from_position,
256                    failure.from_t, failure.seed_uv, failure.error
257                ));
258                if budget.largest() * 2 > CHAIN_CARVE_MAX_PER_SEGMENT {
259                    break Err(failure.error);
260                }
261                budget = budget.doubled();
262            }
263            other => break other,
264        }
265    };
266    if let Err(error) = &first {
267        crate::blend::carve::carve_trace(format_args!(
268            "blend chain: the march at {} refused: {error}", budget.describe()
269        ));
270    }
271    match first {
272        Ok(Rung::Built(built)) => Ok(built),
273        Ok(Rung::TooCoarse { .. }) | Ok(Rung::StationFailed(_)) => {
274            Err(KernelRefusal::internal(KernelStage::Refine, "chain_ladder_rung", "blend: the chain ladder left an accuracy rung unhandled"))
275        }
276        // A fold that reaches a RAIL is not a lens: every section of that wall
277        // is singular somewhere between its contacts, so there is nothing to cut
278        // away and keep, and the ladder below would only re-march it to its top
279        // before the carve refused the same band. Its refusal names the face it
280        // cannot cross, and it is terminal as it stands (`blend/fold.rs`).
281        Err(error) if crate::blend::fold::is_rail_fold(&error) => Err(error),
282        // A CHAMFER section is a straight line and has no ball centre to read
283        // the fold locus against, so its fold stays the refusal it was.
284        Err(error) if crate::blend::is_wall_fold(&error) && !chamfer => {
285            carve_ladder(solid, segments, radius, name, bar)
286                // A CARVE THAT FAILS IS STILL A FOLD, and the fold is terminal.
287                // If the carve's own refusal did not read as one, the ladder
288                // below would answer a proven-folding centre curve with the
289                // cutter's unrelated complaint — which is the exact failure
290                // `fold.rs` was made terminal to stop. So every exit from this
291                // branch carries the fold's own prefix, with the carve's reason
292                // inside it.
293                .map_err(|carve_error| {
294                    if crate::blend::is_wall_fold(&carve_error) {
295                        carve_error
296                    } else {
297                        crate::blend::fold::wall_fold_refusal(
298                            Vec::new(),
299                            format!(
300                                "{} {} fits this edge: the wall folds, and the fold band could not \
301                                 be carved — {carve_error}",
302                                crate::blend::WALL_FOLDS,
303                                radius.abs()
304                            ),
305                        )
306                    }
307                })
308        }
309        Err(error) => Err(error),
310    }
311}
312
313/// Re-march a folding chain up a ladder of station budgets and carve at the
314/// first rung whose wall IS the rolling ball's to `bar`.
315///
316/// The budget is not chosen: a larger or a tighter fillet needs a different
317/// rung, and a number picked once from one fixture would be a resolution
318/// request in disguise. Each rung is measured by the rolling ball itself — its
319/// contacts re-solved halfway between stations, where an interpolant is worst —
320/// against the fitted rails, and the ladder stops at the first rung inside
321/// `intersection_fit`. A wall that still misses at the top is refused by name.
322fn carve_ladder(
323    solid: &BrepSolid,
324    segments: &[ChainSegment<'_>],
325    radius: f64,
326    name: Option<&str>,
327    bar: f64,
328) -> Result<BrepSolid, KernelRefusal> {
329    let mut budget = ChainBudget::for_chain(segments, CHAIN_PER_SEGMENT)?;
330    loop {
331        let per_segment = budget.rung;
332        match build_closed_smooth_chain(
333            solid,
334            segments,
335            radius,
336            false,
337            name,
338            &budget,
339            FoldPolicy::Carve,
340            bar,
341        )? {
342            Rung::Built(built) => {
343                crate::blend::carve::carve_trace(format_args!(
344                    "blend carve: {} is the first rung inside \
345                     {bar:.3e}",
346                    budget.describe()
347                ));
348                return Ok(built);
349            }
350            // The carve keeps its own contract: a station failure refuses as
351            // it did before the plain ladder learned to climb on it.
352            Rung::StationFailed(failure) => return Err(failure.error),
353            Rung::TooCoarse { deviation } => {
354                crate::blend::carve::carve_trace(format_args!(
355                    "blend carve: {} leaves the rails \
356                     {deviation:.6e} from the rolling ball, against {bar:.3e}",
357                    budget.describe()
358                ));
359                if budget.largest() * 2 > CHAIN_CARVE_MAX_PER_SEGMENT {
360                    return Err(crate::blend::fold::wall_fold_refusal(
361                        Vec::new(),
362                        format!(
363                            "{} {} fits this edge: the wall folds, and at {per_segment} stations \
364                             per segment — the top of the carve's ladder — its fitted rails are \
365                             still {deviation:.6e} from the rolling ball's own contacts against a \
366                             bar of {bar:.3e}, so there is no wall accurate enough to carve",
367                            crate::blend::WALL_FOLDS,
368                            radius.abs()
369                        ),
370                    ));
371                }
372                budget = budget.doubled();
373            }
374        }
375    }
376}
377
378/// One rung of a chain build: the finished solid, or — for a carve — the
379/// measured reason this budget is not enough.
380enum Rung {
381    Built(BrepSolid),
382    TooCoarse { deviation: f64 },
383    /// A station's continuation Newton found no root at this budget
384    /// (`march::ChainStationFailure`); a finer rung is a nearer seed.
385    StationFailed(ChainStationFailure),
386}
387
388/// The worst distance from the rolling ball's halfway contacts to the fitted
389/// rails `cr` and `cs`, each found by a closest-point Newton from the
390/// parameter the chord map assigns the midpoint.
391///
392/// THE NEWTON STAYS WHERE IT WAS SEEDED. The foot of a halfway contact is
393/// inside the station interval it was solved in, so every step is bounded by
394/// that interval's own length in the fit's parameter, clamped to the rail's
395/// domain, and taken only when it brings the rail closer. A plain Newton stops
396/// wherever the distance is stationary along the curve, which a rail that
397/// turns sharply at a chain junction offers well away from the foot. Measured
398/// on a 20 cube rounded at r = 3 with a bottom edge at r = 4 (the wall folds at
399/// the arcs' rails, so this is the ladder's reading with that refusal set
400/// aside): unguarded, every rung read the distance from the halfway contact
401/// beside a junction to the junction itself — 2.927709e-1 at 24 stations per
402/// segment against a half interval of 0.291667, halving per rung to
403/// 9.128591e-3 at 768 — and on two of the eight mirror-image edges a foot
404/// across the chain, 1.299089e1. Guarded, all eight read 6.961233e-3 at 24
405/// and 6.920006e-6 to 6.920008e-6 at 768, falling fourfold per rung as a
406/// cubic's does. A far foot can read short as well as long, which would pass
407/// a rung that misses.
408///
409/// Also returns every interval over `bar`, named by its left sample's index:
410/// what local refinement inserts a station into.
411fn rail_deviation(
412    samples: &[ChainSample],
413    cr: &NurbsCurve,
414    cs: &NurbsCurve,
415    fit_low: f64,
416    fit_range: f64,
417    budget: &ChainBudget,
418    bar: f64,
419) -> Result<(f64, Vec<usize>), KernelRefusal> {
420    // Each interval runs from a station to the NEXT one along its segment — the
421    // next grid position, a station local refinement inserted, or the
422    // segment's own end station — in the order the stations take.
423    let mut order: Vec<usize> = (0..samples.len())
424        .filter(|&index| (0..=budget.count(samples[index].segment) as isize).contains(&samples[index].position))
425        .collect();
426    order.sort_by(|&a, &b| {
427        let (a, b) = (&samples[a], &samples[b]);
428        a.segment.cmp(&b.segment).then(a.position.cmp(&b.position)).then(a.offset.total_cmp(&b.offset))
429    });
430    let next_of: std::collections::HashMap<usize, f64> = order
431        .windows(2)
432        .filter(|pair| samples[pair[0]].segment == samples[pair[1]].segment)
433        .map(|pair| (pair[0], samples[pair[1]].parameter))
434        .collect();
435    let mut failing = Vec::new();
436    let closest = closest_on_rail;
437    let mut worst: f64 = 0.0;
438    // Each segment's worst halfway contact, for `BREP_BLEND_CARVE_TRACE=1`: where
439    // along the chain a rung misses is what tells a junction from a bend.
440    let mut per_segment_worst: Vec<(f64, isize, f64, f64, Vec3, Vec3)> = Vec::new();
441    for &index in &order {
442        let sample = &samples[index];
443        let Some([p1, p2]) = sample.midpoint else {
444            continue;
445        };
446        let Some(&next) = next_of.get(&index) else {
447            continue;
448        };
449        let global = (0.5 * (sample.parameter + next)).rem_euclid(1.0);
450        let seed = (global - fit_low) / fit_range;
451        let reach = (next - sample.parameter).abs() / fit_range;
452        let (miss1, miss2) = (closest(cr, p1, seed, reach)?, closest(cs, p2, seed, reach)?);
453        worst = worst.max(miss1).max(miss2);
454        if miss1.max(miss2) > bar {
455            failing.push(index);
456        }
457        if per_segment_worst.len() <= sample.segment {
458            per_segment_worst.resize(
459                sample.segment + 1,
460                (0.0, 0, 0.0, 0.0, Vec3::default(), Vec3::default()),
461            );
462        }
463        if miss1.max(miss2) >= per_segment_worst[sample.segment].0 {
464            per_segment_worst[sample.segment] = (miss1.max(miss2), sample.position, miss1, miss2, p1, p2);
465        }
466    }
467    for (segment, (miss, position, miss1, miss2, p1, p2)) in per_segment_worst.iter().enumerate() {
468        crate::blend::carve::carve_trace(format_args!(
469            "  rail miss segment {segment}: worst {miss:.6e} after station {position} \
470             (cr {miss1:.6e} at ({:.6}, {:.6}, {:.6}), cs {miss2:.6e} at ({:.6}, {:.6}, {:.6}))",
471            p1.x, p1.y, p1.z, p2.x, p2.y, p2.z
472        ));
473    }
474    Ok((worst, failing))
475}
476
477/// Most dense INTERVALS (4097 samples) a closed chain's support-piece track may
478/// hold: the existing ceiling its uniform track is built under, and the one
479/// local refinement of that track stays within — not a sample count.
480const PIECE_TRACK_CEILING: usize = 4096;
481
482/// Distance from `target` to a fitted rail, by a Newton that stays where it is
483/// seeded (`rail_deviation`'s reading, shared with the dense diagnostic).
484fn closest_on_rail(curve: &NurbsCurve, target: Vec3, seed: f64, reach: f64) -> Result<f64, KernelRefusal> {
485    let [low, high] = curve.domain().or_refuse(KernelStage::Refine, "domain")?;
486    let distance = |u: f64| -> Result<f64, KernelRefusal> {
487        Ok(curve.derivatives_extended(u, 0).or_refuse(KernelStage::Refine, "derivatives_extended")?[0].sub(target).length())
488    };
489    let mut u = seed.clamp(low, high);
490    let mut best = distance(u)?;
491    for _ in 0..40 {
492        let derivatives = curve.derivatives_extended(u, 2).or_refuse(KernelStage::Refine, "derivatives_extended")?;
493        let offset = derivatives[0].sub(target);
494        let slope = offset.dot(derivatives[1]);
495        let curvature = derivatives[1].dot(derivatives[1]) + offset.dot(derivatives[2]);
496        // Off a minimum's basin the Newton step points the wrong way; walk
497        // downhill by the bound instead and let the test below shorten it.
498        let mut step = if curvature > 0.0 {
499            (slope / curvature).clamp(-reach, reach)
500        } else {
501            reach.copysign(slope)
502        };
503        let mut trial = (u - step).clamp(low, high);
504        let mut trial_distance = distance(trial)?;
505        for _ in 0..30 {
506            if trial_distance < best {
507                break;
508            }
509            step *= 0.5;
510            trial = (u - step).clamp(low, high);
511            trial_distance = distance(trial)?;
512        }
513        if !(trial_distance < best) {
514            break;
515        }
516        let moved = (trial - u).abs();
517        u = trial;
518        best = trial_distance;
519        if moved < 1e-15 {
520            break;
521        }
522    }
523    Ok(best)
524}
525
526/// Points of every interval the LOCAL-INSERTION path judges a chain rung's
527/// rails on: j/16, j = 1..15 (the halfway point among them).
528pub(super) const LOCAL_RAIL_STENCIL: usize = 16;
529
530/// One dense contact pair's miss: each rail's raw reading is checked BEFORE
531/// they are joined (`f64::max` alone drops a NaN on either side); `None` when
532/// either is unreadable.
533fn rail_pair_miss(cr_miss: f64, cs_miss: f64) -> Option<f64> {
534    (cr_miss.is_finite() && cs_miss.is_finite()).then(|| cr_miss.max(cs_miss))
535}
536
537/// Each solved dense contact's distance to the fitted rails, read with the same
538/// closest-point Newton and seeds `rail_deviation` uses, keyed by the
539/// interval's left sample: the readable misses, and separately every interval
540/// with an unreadable reading on either rail (a non-finite distance, or a foot
541/// that could not be evaluated).
542#[allow(clippy::too_many_arguments)]
543pub(super) fn dense_rail_misses(
544    contacts: &crate::blend::chain::DenseContacts,
545    samples: &[ChainSample],
546    cr: &NurbsCurve,
547    cs: &NurbsCurve,
548    low: f64,
549    range: f64,
550    wrap: bool,
551) -> Result<(Vec<(usize, f64)>, Vec<usize>), KernelRefusal> {
552    let mut readable = Vec::with_capacity(contacts.solved.len());
553    let mut unreadable = Vec::new();
554    for (left, next, fraction, p1, p2) in &contacts.solved {
555        let (from, to) = (samples[*left].parameter, samples[*next].parameter);
556        // A closed chain's chord parameter wraps with period 1; an OPEN
557        // chain's runs straight through its ends (its last station's may
558        // overshoot 1, and wrapping it would read the far end).
559        let global = from + fraction * (to - from);
560        let global = if wrap { global.rem_euclid(1.0) } else { global };
561        let seed = (global - low) / range;
562        let reach = (to - from).abs() / range;
563        // A closest foot that cannot be evaluated is an unreadable reading of
564        // this interval, not a refusal of the rung.
565        let cr_miss = closest_on_rail(cr, *p1, seed, reach).unwrap_or(f64::NAN);
566        let cs_miss = closest_on_rail(cs, *p2, seed, reach).unwrap_or(f64::NAN);
567        match rail_pair_miss(cr_miss, cs_miss) {
568            Some(miss) => readable.push((*left, miss)),
569            None => unreadable.push(*left),
570        }
571    }
572    Ok((readable, unreadable))
573}
574
575/// The local-insertion path's rail verdict: the halfway reading joined with the
576/// dense stencil's, FAIL-CLOSED — an interval with a dense contact over `bar`,
577/// an unreadable distance, or a contact that was not read (unsolved, or
578/// unreadable on either rail), fails; the worst is the largest readable
579/// distance. Intervals are named by their left sample.
580pub(super) fn combine_rail_verdict(
581    halfway: (f64, Vec<usize>),
582    dense: &[(usize, f64)],
583    unsolved: &[usize],
584    bar: f64,
585) -> (f64, Vec<usize>) {
586    let (mut worst, mut failing) = halfway;
587    for &(left, miss) in dense {
588        if !miss.is_finite() {
589            failing.push(left);
590            continue;
591        }
592        if miss > bar {
593            failing.push(left);
594        }
595        worst = worst.max(miss);
596    }
597    failing.extend_from_slice(unsolved);
598    failing.sort_unstable();
599    failing.dedup();
600    (worst, failing)
601}
602
603/// Whether the construction request runs on an accepted plain fillet rung:
604/// always, outside the controls that turn it off on their own thread
605/// (`MODEL_REQUEST_OFF`) or that turn the wall verdict off.
606/// The construction request's round limit: the existing local-refinement
607/// budget, outside the control that sets it to zero on its own thread to
608/// force the typed unmet report.
609pub(super) fn model_rounds_limit() -> usize {
610    REFINE_ROUNDS
611}
612
613pub(super) fn model_request_on() -> bool {
614    wall_verdict_on()
615}
616
617
618/// Whether the fillet wall verdict runs: always, outside the mutation control
619/// that turns it off on its own thread (`WALL_VERDICT_OFF`).
620pub(super) fn wall_verdict_on() -> bool {
621    true
622}
623
624
625
626
627
628/// Section-angle probes across a fillet wall, as fractions of its `v` span:
629/// k/8 for k = 1..7, between (never on) the rails the rail verdict reads.
630pub(super) const WALL_SECTION_PROBES: usize = 8;
631
632/// How far `p` stands from the rolling ball's sweep, read at p's OWN
633/// characteristic: the stationary point of |p − C(s)|, C the cubic through
634/// four consecutive exact ball centres of `nodes` around node `j` (`nodes`
635/// uniform in fraction across one interval, its ends the interval's own
636/// stations), searched within one node of `j`, minus `radius`.
637///
638/// The cubic is an APPROXIMATION of the centre curve, not the exact sweep:
639/// its error grows with the node spacing and the spine's higher derivatives,
640/// and [`wall_probe_reading`] estimates it by a second, overlapping stencil.
641///
642/// A STATIONARY point, not the nearest centre: past a fold's focal point an
643/// envelope point is a MAXIMUM of its distance to the spine, and a global
644/// nearest centre would read a carved wall's lens as inside another ball. A
645/// centre at the probe's own chain parameter is not used either: the wall's
646/// `u` and the stencil's edge parameter name the same section only to first
647/// order, and a spine offset ℓ adds √(ρ² + ℓ²) − ρ to the reading.
648///
649/// `None` — unreadable, which fails its interval — for a non-finite input,
650/// fewer than four nodes, a probe at an interval end, or a solve that leaves
651/// its window or does not converge.
652pub(in crate::blend) fn sweep_residual(p: Vec3, nodes: &[Vec3], j: usize, radius: f64) -> Option<f64> {
653    let j = locate_foot(p, nodes, j, radius)?;
654    sweep_residual_on(p, nodes, (j - 1).min(nodes.len() - 4), j, radius)
655}
656
657/// How many nodes the stencil may follow a probe's foot away from the probe's
658/// own node.
659const FOOT_TRAVEL: usize = 4;
660
661/// The node whose stencil holds `p`'s stationary point, starting from node
662/// `j`: the probe's wall `u` and the centres' edge parameter name the same
663/// section only to first order, so the foot can sit more than the one-node
664/// window from `j` (the 15- and 20-degree rows' feet stood up to 1.5 nodes,
665/// one 3.3, from it on f4ed17782's build, every decline there a window exit).
666/// A solve that leaves its window moves the stencil to the node nearest
667/// where it left, at most [`FOOT_TRAVEL`] nodes from `j` and never back to a
668/// node already tried; `None` when the foot is not found that way (or a
669/// solve is unreadable), which fails its interval as before.
670fn locate_foot(p: Vec3, nodes: &[Vec3], j: usize, radius: f64) -> Option<usize> {
671    locate_foot_read(p, nodes, j, radius).map(|(node, _)| node)
672}
673
674/// Keep the converged stencil's distance with its node. The caller can use
675/// this exact reading instead of solving the same stencil a second time.
676fn locate_foot_read(p: Vec3, nodes: &[Vec3], j: usize, radius: f64) -> Option<(usize, f64)> {
677    locate_foot_read_prepared(p, nodes, j, radius, None)
678}
679
680fn locate_foot_read_prepared(p: Vec3, nodes: &[Vec3], j: usize, radius: f64, prepared: Option<&StencilPreparation>) -> Option<(usize, f64)> {
681    let n = nodes.len();
682    if n < 4 || j == 0 || j + 1 >= n {
683        return None;
684    }
685    let mut tried: Vec<usize> = Vec::new();
686    let mut at = j;
687    loop {
688        let base = (at - 1).min(n - 4);
689        match stencil_solve_prepared(p, nodes, base, at, radius, prepared) {
690            Solve::Read(distance) => return Some((at, distance)),
691            Solve::Unread => return None,
692            Solve::Window(s) => {
693                tried.push(at);
694                let next = (base as f64 + s).round();
695                let next = if next.is_finite() { next.clamp(1.0, (n - 2) as f64) as usize } else { at };
696                if tried.contains(&next) || next.abs_diff(j) > FOOT_TRAVEL {
697                    wall_decline(|| format!(
698                        "foot search from node {j} of {n}: tried {tried:?}, next {next} ({})",
699                        if tried.contains(&next) { "already tried" } else { "beyond FOOT_TRAVEL" }
700                    ));
701                    return None;
702                }
703                at = next;
704            }
705        }
706    }
707}
708
709/// [`sweep_residual`] on the cubic through `nodes[base..base + 4]`, the
710/// stationary point searched from node `j` within one node of it.
711fn sweep_residual_on(p: Vec3, nodes: &[Vec3], base: usize, j: usize, radius: f64) -> Option<f64> {
712    match stencil_solve(p, nodes, base, j, radius) {
713        Solve::Read(distance) => Some(distance),
714        Solve::Window(_) | Solve::Unread => None,
715    }
716}
717
718/// One stencil's solve: the distance read, or where (in node units from
719/// `base`) the search left its one-node window, or unreadable.
720enum Solve {
721    Read(f64),
722    Window(f64),
723    Unread,
724}
725
726thread_local! {
727    /// When Some, every decline inside the declared wall reader appends WHY
728    /// (innermost first): set and taken by `declared_wall_probe_why` only, so
729    /// every other caller pays one thread-local read per decline and the
730    /// reader's arithmetic and outcomes are the same either way.
731    static WALL_DECLINES: std::cell::RefCell<Option<Vec<String>>> = const { std::cell::RefCell::new(None) };
732}
733
734fn wall_decline(reason: impl FnOnce() -> String) {
735    WALL_DECLINES.with(|log| {
736        if let Some(lines) = log.borrow_mut().as_mut() {
737            lines.push(reason());
738        }
739    });
740}
741
742/// [`declared_wall_probe`] with, when it declines, the chain of declines
743/// that made it (innermost first): which tier (fine / coarse), which
744/// stencil, and inside the solver the exit (unreadable input, a vanishing
745/// slope, a window exit with the overshot `s`, its offset from the window and
746/// whether the solver's OWN stationarity predicate already held at the
747/// overshot point, a stall, a non-finite distance) or the foot search's
748/// tried nodes. Exactly the same arithmetic as `declared_wall_probe`.
749pub(in crate::blend) fn declared_wall_probe_why(p: Vec3, fine: &[Vec3], jf: usize, radius: f64) -> Result<(f64, f64), String> {
750    let previous = WALL_DECLINES.with(|log| log.replace(Some(Vec::new())));
751    let read = declared_wall_probe(p, fine, jf, radius);
752    let lines = WALL_DECLINES.with(|log| log.replace(previous)).unwrap_or_default();
753    read.ok_or_else(|| if lines.is_empty() { "declined (no inner cause recorded)".to_string() } else { lines.join(" <- ") })
754}
755
756fn stencil_solve(p: Vec3, nodes: &[Vec3], base: usize, j: usize, radius: f64) -> Solve {
757    stencil_solve_prepared(p, nodes, base, j, radius, None)
758}
759
760fn stencil_solve_prepared(p: Vec3, nodes: &[Vec3], base: usize, j: usize, radius: f64, prepared: Option<&StencilPreparation>) -> Solve {
761    let finite = |v: Vec3| v.x.is_finite() && v.y.is_finite() && v.z.is_finite();
762    let n = nodes.len();
763    if n < 4 || j == 0 || j + 1 >= n || base + 4 > n || !radius.is_finite() || !finite(p) || !prepared.map_or_else(|| wall_nodes_finite(nodes), |cache| cache.finite) {
764        wall_decline(|| format!("stencil {n} nodes base {base} j {j}: unreadable input"));
765        return Solve::Unread;
766    }
767    let [p0, d1, d2, d3] = prepared.map_or_else(|| wall_coefficients(nodes, base), |cache| cache.coefficients(nodes, base));
768    let curve = |s: f64| -> (Vec3, Vec3, Vec3) {
769        let c = p0.add(d1.scale(s)).add(d2.scale(s * (s - 1.0))).add(d3.scale(s * (s - 1.0) * (s - 2.0)));
770        let c1 = d1.add(d2.scale(2.0 * s - 1.0)).add(d3.scale(3.0 * s * s - 6.0 * s + 2.0));
771        let c2 = d2.scale(2.0).add(d3.scale(6.0 * s - 6.0));
772        (c, c1, c2)
773    };
774    let start = j as f64 - base as f64;
775    let mut s = start;
776    let mut converged = false;
777    for _ in 0..50 {
778        let (c, c1, c2) = curve(s);
779        let w = c.sub(p);
780        let (g, slope) = (w.dot(c1), c1.dot(c1) + w.dot(c2));
781        if !(g.is_finite() && slope.is_finite()) || slope == 0.0 {
782            wall_decline(|| format!("stencil {n} nodes base {base} j {j}: slope {slope:e} g {g:e} at s {s}"));
783            return Solve::Unread;
784        }
785        // Converged at THIS s once the stationarity condition holds to its own
786        // rounding: g = (C − p)·C′ carries an absolute error of a few
787        // ε·(|p| + |C|)·|C′|, so at a node spacing of 1e-4 the Newton step
788        // cannot fall below ~1e-11 node units and the relative step test
789        // alone never fires — a line-for-line model of this reader read 87 %
790        // to 97 % of the NX and skew walls' fine-stencil probes unreadable
791        // that way (its Newton run past the stop read |g| / (ε·(|p| + |C|)·|C′|)
792        // at most 0.91 there; 16 keeps a margin). Tested before the step, so
793        // the point read is the point whose stationarity was checked; its
794        // distance is off the true stationary one by about
795        // ½·|g|·δs/|C − p| ≤ 8ε·(|p| + |C|)·|C′|/ρ for δs within the
796        // one-node window, far below any bar.
797        let stationary = g.abs() <= 16.0 * f64::EPSILON * (p.length() + c.length()) * c1.length();
798        if stationary {
799            converged = true;
800            break;
801        }
802        let step = g / slope;
803        let before = s;
804        s -= step;
805        if !((s - start).abs() <= 1.0) {
806            wall_decline(|| {
807                // The solver's own stationarity predicate, read at the
808                // overshot point (diagnostic only: nothing decides on it).
809                let (c, c1, _) = curve(s);
810                let g_after = c.sub(p).dot(c1);
811                let stationary_after = g_after.abs() <= 16.0 * f64::EPSILON * (p.length() + c.length()) * c1.length();
812                format!(
813                    "stencil {n} nodes base {base} j {j}: window exit s {before} -> {s} (start {start}, |s - start| - 1 = {:e}), \
814                     g {g:e} -> {g_after:e}, stationary at the overshot s: {stationary_after}",
815                    (s - start).abs() - 1.0
816                )
817            });
818            return if s.is_finite() { Solve::Window(s) } else { Solve::Unread };
819        }
820        if step.abs() <= 1e-14 * (1.0 + s.abs()) {
821            converged = true;
822            break;
823        }
824    }
825    if !converged {
826        wall_decline(|| format!("stencil {n} nodes base {base} j {j}: no convergence in 50 steps, s {s}"));
827        return Solve::Unread;
828    }
829    let distance = curve(s).0.sub(p).length() - radius;
830    if distance.is_finite() {
831        Solve::Read(distance)
832    } else {
833        wall_decline(|| format!("stencil {n} nodes base {base} j {j}: non-finite distance at s {s}"));
834        Solve::Unread
835    }
836}
837
838/// One wall probe read on TWO overlapping four-centre stencils around node
839/// `j` (the one [`sweep_residual`] uses and its neighbour): `(worst, spread)`
840/// — the larger magnitude, which keeps ONE stencil's under-read from passing,
841/// and how far apart the two readings are. Both are PROXY readings: when both
842/// stencils under-read alike (a symmetric node set about a foot near the
843/// focal distance) the larger still under-reads and the spread reads zero,
844/// so neither number is an error bound. `None` when either stencil is
845/// unreadable.
846pub(in crate::blend) fn wall_probe_reading(p: Vec3, nodes: &[Vec3], j: usize, radius: f64) -> Option<(f64, f64)> {
847    let n = nodes.len();
848    if n < 5 || j == 0 || j + 1 >= n {
849        return None;
850    }
851    let j = locate_foot(p, nodes, j, radius)?;
852    let first = (j - 1).min(n - 4);
853    let second = if j >= 2 { (j - 2).min(n - 4) } else { j.min(n - 4) };
854    let second = if second == first { if first >= 1 { first - 1 } else { first + 1 } } else { second };
855    let a = sweep_residual_on(p, nodes, first, j, radius)?;
856    let b = sweep_residual_on(p, nodes, second, j, radius)?;
857    Some((a.abs().max(b.abs()), (a - b).abs()))
858}
859
860/// The DECLARED wall reading at one probe, the one rule both wall tiers
861/// decide on: `(reading, uncertainty)`, declared = reading + uncertainty.
862/// `fine` are the interval's centres at twice the probe density, `jf` the
863/// probe's (even) node among them; the coarse set is every other node. The
864/// reading is the larger fine two-stencil magnitude; the uncertainty is the
865/// larger of the fine stencils' spread and how far the fine reading moved
866/// from the coarse one — the proxy's error estimated from its own change with
867/// node density, which two stencils agreeing at one density cannot show (two
868/// symmetric stencils that both under-read read a spread of zero). An
869/// ESTIMATE, not a bound: a proxy whose error does not shrink with density
870/// would still hide. `None` when any reading is unreadable.
871
872fn wall_solve_reuse_on() -> bool {
873    std::env::var("BREP_WALL_SOLVE_REUSE").map_or(true, |value| value != "0" && value != "off")
874}
875
876pub(in crate::blend) fn declared_wall_probe(p: Vec3, fine: &[Vec3], jf: usize, radius: f64) -> Option<(f64, f64)> {
877    declared_wall_probe_with(p, fine, jf, radius, wall_solve_reuse_on())
878}
879
880fn declared_wall_probe_with(p: Vec3, fine: &[Vec3], jf: usize, radius: f64, reuse: bool) -> Option<(f64, f64)> {
881    declared_wall_probe_prepared(p, fine, jf, radius, reuse, None)
882}
883
884fn declared_wall_probe_prepared(p: Vec3, fine: &[Vec3], jf: usize, radius: f64, reuse: bool, prepared: Option<&PreparedWall<'_>>) -> Option<(f64, f64)> {
885    if jf % 2 != 0 || fine.len() < 9 || fine.len() % 2 == 0 {
886        wall_decline(|| format!("declared reader: node {jf} of {} not an even node of an odd fine set", fine.len()));
887        return None;
888    }
889    let (reading, spread, at_fine) = if reuse {
890        let read = (|| {
891            let n = fine.len();
892            let (j, a) = locate_foot_read_prepared(p, fine, jf, radius, prepared.map(|c| &c.fine_stencils))?;
893            let first = (j - 1).min(n - 4);
894            let second = if j >= 2 { (j - 2).min(n - 4) } else { j.min(n - 4) };
895            let second = if second == first { if first >= 1 { first - 1 } else { first + 1 } } else { second };
896            let b = match stencil_solve_prepared(p, fine, second, j, radius, prepared.map(|c| &c.fine_stencils)) {
897                Solve::Read(distance) => distance,
898                _ => return None,
899            };
900            Some((a.abs().max(b.abs()), (a - b).abs(), a))
901        })();
902        let Some(read) = read else {
903            wall_decline(|| format!("FINE two-stencil reading at node {jf} of {}", fine.len()));
904            return None;
905        };
906        read
907    } else {
908        let Some((reading, spread)) = wall_probe_reading(p, fine, jf, radius) else {
909            wall_decline(|| format!("FINE two-stencil reading at node {jf} of {}", fine.len()));
910            return None;
911        };
912        let Some(at_fine) = sweep_residual(p, fine, jf, radius) else {
913            wall_decline(|| format!("FINE sweep at node {jf} of {}", fine.len()));
914            return None;
915        };
916        (reading, spread, at_fine)
917    };
918    let coarse_owned;
919    let coarse = match prepared {
920        Some(cache) => cache.coarse.as_slice(),
921        None => { coarse_owned = wall_coarse(fine); &coarse_owned }
922    };
923    let at_coarse = if reuse {
924        locate_foot_read_prepared(p, coarse, jf / 2, radius, prepared.map(|c| &c.coarse_stencils)).map(|(_, distance)| distance)
925    } else {
926        sweep_residual(p, coarse, jf / 2, radius)
927    };
928    let Some(at_coarse) = at_coarse else {
929        wall_decline(|| format!("COARSE sweep at node {} of {}", jf / 2, coarse.len()));
930        return None;
931    };
932    let uncertainty = spread.max((at_fine - at_coarse).abs());
933    uncertainty.is_finite().then_some((reading, uncertainty))
934}
935
936/// The even fine node whose position along the CENTRE PATH (cumulative
937/// chord of `fine`) is nearest `fraction` of it, inside 2..=n-3; `None` when
938/// the path has no finite positive length.
939///
940/// A wall probe sits at a fraction of its interval in the WALL's parameter,
941/// which is chord-like; the exact centres sit at uniform fractions of the
942/// edge's NATIVE parameter. On a regular edge the two agree to first order and
943/// the probe's foot lies within a node or two of its own node. At a
944/// STATIONARY section (the edge's derivative vanishes there) the centre moves
945/// like t^2 in the native parameter, so a probe at chord fraction f has its foot
946/// near native fraction sqrt(f): the probe at node 2 of 32 (f = 1/16) has its
947/// foot near node 8, past `FOOT_TRAVEL`, and the reader declines though the
948/// centre locus is perfectly regular there (zero_tangent_sites' stationary
949/// stadium and tombstone rims refused that way with every reading under the
950/// bar).
951pub(in crate::blend) fn centre_path_node(fine: &[Vec3], fraction: f64) -> Option<usize> {
952    let n = fine.len();
953    if n < 7 || !(fraction.is_finite()) {
954        return None;
955    }
956    let mut cumulative = Vec::with_capacity(n);
957    cumulative.push(0.0_f64);
958    for pair in fine.windows(2) {
959        let step = pair[1].sub(pair[0]).length();
960        cumulative.push(cumulative.last().copied().unwrap_or(0.0) + step);
961    }
962    let total = *cumulative.last()?;
963    if !(total.is_finite() && total > 0.0) {
964        return None;
965    }
966    let target = fraction.clamp(0.0, 1.0) * total;
967    let nearest = (0..n).min_by(|&a, &b| (cumulative[a] - target).abs().total_cmp(&(cumulative[b] - target).abs()))?;
968    let even = 2 * ((nearest as f64 / 2.0).round() as usize);
969    Some(even.clamp(2, (n - 3) & !1))
970}
971
972/// [`declared_wall_probe_why`] from the probe's own node `jf`; when that
973/// declines, once more from the node the probe's `fraction` reaches along the
974/// centre path ([`centre_path_node`]), when that is another node. The same
975/// reader, bars and fail-closed rule decide; a probe it still cannot read is
976/// unread (both declines named). A probe the first read reads is read exactly
977/// as before.
978pub(in crate::blend) fn declared_wall_probe_located_why(p: Vec3, fine: &[Vec3], jf: usize, fraction: f64, radius: f64) -> Result<(f64, f64), String> {
979    let first = match declared_wall_probe_why(p, fine, jf, radius) {
980        Ok(read) => return Ok(read),
981        Err(why) => why,
982    };
983    match centre_path_node(fine, fraction) {
984        Some(located) if located != jf => declared_wall_probe_why(p, fine, located, radius)
985            .map_err(|again| format!("{first} -- and from centre-path node {located}: {again}")),
986        _ => Err(first),
987    }
988}
989
990/// [`declared_wall_probe_located_why`] as an `Option`.
991pub(in crate::blend) fn declared_wall_probe_located(p: Vec3, fine: &[Vec3], jf: usize, fraction: f64, radius: f64) -> Option<(f64, f64)> {
992    match declared_wall_probe(p, fine, jf, radius) {
993        Some(read) => Some(read),
994        None => match centre_path_node(fine, fraction) {
995            Some(located) if located != jf => declared_wall_probe(p, fine, located, radius),
996            _ => None,
997        },
998    }
999}
1000
1001/// [`declared_wall_probe_why`] on centres at NON-UNIFORM native abscissae
1002/// `ts` (one per node, strictly increasing): the same two-tier rule, the
1003/// same two overlapping four-centre stencils, the same foot search
1004/// (`FOOT_TRAVEL`), the same 16·ε stationarity predicate and the same
1005/// window (one node either side of the start node), with each cubic the
1006/// Newton interpolant on the nodes' OWN abscissae instead of the implicit
1007/// 0, 1, 2, 3. For a junction of two intervals of unequal native width (local
1008/// insertions make them so): the uniform reader would place the two halves'
1009/// nodes at one spacing and distort the centre locus. The coarse tier is
1010/// every other node at every other abscissa. A caller whose abscissae are
1011/// uniform uses the uniform reader, bit for bit.
1012pub(in crate::blend) fn declared_wall_probe_at_why(p: Vec3, fine: &[Vec3], ts: &[f64], jf: usize, radius: f64) -> Result<(f64, f64), String> {
1013    if jf % 2 != 0
1014        || fine.len() < 9
1015        || fine.len() % 2 == 0
1016        || ts.len() != fine.len()
1017        || !ts.iter().all(|t| t.is_finite())
1018        || !ts.windows(2).all(|pair| pair[1] > pair[0])
1019    {
1020        return Err(format!("native-abscissa reader: node {jf} of {} with {} finite, strictly increasing abscissae required", fine.len(), ts.len()));
1021    }
1022    let previous = WALL_DECLINES.with(|log| log.replace(Some(Vec::new())));
1023    let read = (|| -> Option<(f64, f64)> {
1024        let (reading, spread) = wall_probe_reading_at(p, fine, ts, jf, radius).or_else(|| {
1025            wall_decline(|| format!("FINE two-stencil reading (native abscissae) at node {jf} of {}", fine.len()));
1026            None
1027        })?;
1028        let at_fine = sweep_residual_at(p, fine, ts, jf, radius).or_else(|| {
1029            wall_decline(|| format!("FINE sweep (native abscissae) at node {jf} of {}", fine.len()));
1030            None
1031        })?;
1032        let coarse: Vec<Vec3> = fine.iter().step_by(2).copied().collect();
1033        let coarse_ts: Vec<f64> = ts.iter().step_by(2).copied().collect();
1034        let at_coarse = sweep_residual_at(p, &coarse, &coarse_ts, jf / 2, radius).or_else(|| {
1035            wall_decline(|| format!("COARSE sweep (native abscissae) at node {} of {}", jf / 2, coarse.len()));
1036            None
1037        })?;
1038        let uncertainty = spread.max((at_fine - at_coarse).abs());
1039        uncertainty.is_finite().then_some((reading, uncertainty))
1040    })();
1041    let lines = WALL_DECLINES.with(|log| log.replace(previous)).unwrap_or_default();
1042    read.ok_or_else(|| if lines.is_empty() { "declined (no inner cause recorded)".to_string() } else { lines.join(" <- ") })
1043}
1044
1045fn wall_probe_reading_at(p: Vec3, nodes: &[Vec3], ts: &[f64], j: usize, radius: f64) -> Option<(f64, f64)> {
1046    let n = nodes.len();
1047    if n < 5 || j == 0 || j + 1 >= n {
1048        return None;
1049    }
1050    let j = locate_foot_at(p, nodes, ts, j, radius)?;
1051    let first = (j - 1).min(n - 4);
1052    let second = if j >= 2 { (j - 2).min(n - 4) } else { j.min(n - 4) };
1053    let second = if second == first { if first >= 1 { first - 1 } else { first + 1 } } else { second };
1054    let a = match stencil_solve_at(p, nodes, ts, first, j, radius) { Solve::Read(distance) => distance, _ => return None };
1055    let b = match stencil_solve_at(p, nodes, ts, second, j, radius) { Solve::Read(distance) => distance, _ => return None };
1056    Some((a.abs().max(b.abs()), (a - b).abs()))
1057}
1058
1059fn sweep_residual_at(p: Vec3, nodes: &[Vec3], ts: &[f64], j: usize, radius: f64) -> Option<f64> {
1060    let j = locate_foot_at(p, nodes, ts, j, radius)?;
1061    match stencil_solve_at(p, nodes, ts, (j - 1).min(nodes.len() - 4), j, radius) {
1062        Solve::Read(distance) => Some(distance),
1063        Solve::Window(_) | Solve::Unread => None,
1064    }
1065}
1066
1067/// [`locate_foot`] on native abscissae: a window exit moves the stencil to
1068/// the node whose abscissa is nearest the exit, within `FOOT_TRAVEL` nodes
1069/// and never back to a node already tried.
1070fn locate_foot_at(p: Vec3, nodes: &[Vec3], ts: &[f64], j: usize, radius: f64) -> Option<usize> {
1071    let n = nodes.len();
1072    if n < 4 || j == 0 || j + 1 >= n {
1073        return None;
1074    }
1075    let mut tried: Vec<usize> = Vec::new();
1076    let mut at = j;
1077    loop {
1078        let base = (at - 1).min(n - 4);
1079        match stencil_solve_at(p, nodes, ts, base, at, radius) {
1080            Solve::Read(_) => return Some(at),
1081            Solve::Unread => return None,
1082            Solve::Window(t) => {
1083                tried.push(at);
1084                let nearest = (0..n).min_by(|&a, &b| (ts[a] - t).abs().total_cmp(&(ts[b] - t).abs())).unwrap_or(at);
1085                let next = if t.is_finite() { nearest.clamp(1, n - 2) } else { at };
1086                if tried.contains(&next) || next.abs_diff(j) > FOOT_TRAVEL {
1087                    wall_decline(|| format!(
1088                        "foot search (native abscissae) from node {j} of {n}: tried {tried:?}, next {next} ({})",
1089                        if tried.contains(&next) { "already tried" } else { "beyond FOOT_TRAVEL" }
1090                    ));
1091                    return None;
1092                }
1093                at = next;
1094            }
1095        }
1096    }
1097}
1098
1099/// [`stencil_solve`] on the four nodes `base..base + 4` at their native
1100/// abscissae: the Newton (divided-difference) cubic in `t`, searched from
1101/// node `j`'s abscissa within one node either side; `Window` carries the
1102/// abscissa it left at.
1103fn stencil_solve_at(p: Vec3, nodes: &[Vec3], ts: &[f64], base: usize, j: usize, radius: f64) -> Solve {
1104    let finite = |v: Vec3| v.x.is_finite() && v.y.is_finite() && v.z.is_finite();
1105    let n = nodes.len();
1106    if n < 4 || j == 0 || j + 1 >= n || base + 4 > n || ts.len() != n || !radius.is_finite() || !finite(p) || !nodes.iter().all(|&c| finite(c)) {
1107        wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: unreadable input"));
1108        return Solve::Unread;
1109    }
1110    let [p0, p1, p2, p3] = [nodes[base], nodes[base + 1], nodes[base + 2], nodes[base + 3]];
1111    let [t0, t1, t2, t3] = [ts[base], ts[base + 1], ts[base + 2], ts[base + 3]];
1112    let d1 = p1.sub(p0).scale(1.0 / (t1 - t0));
1113    let d12 = p2.sub(p1).scale(1.0 / (t2 - t1));
1114    let d23 = p3.sub(p2).scale(1.0 / (t3 - t2));
1115    let d2 = d12.sub(d1).scale(1.0 / (t2 - t0));
1116    let d123 = d23.sub(d12).scale(1.0 / (t3 - t1));
1117    let d3 = d123.sub(d2).scale(1.0 / (t3 - t0));
1118    let curve = |t: f64| -> (Vec3, Vec3, Vec3) {
1119        let (a, b, c) = (t - t0, t - t1, t - t2);
1120        let point = p0.add(d1.scale(a)).add(d2.scale(a * b)).add(d3.scale(a * b * c));
1121        let first = d1.add(d2.scale(a + b)).add(d3.scale(b * c + a * c + a * b));
1122        let second = d2.scale(2.0).add(d3.scale(2.0 * (a + b + c)));
1123        (point, first, second)
1124    };
1125    let (low, high) = (ts[j - 1], ts[j + 1]);
1126    let start = ts[j];
1127    let mut t = start;
1128    let mut converged = false;
1129    for _ in 0..50 {
1130        let (c, c1, c2) = curve(t);
1131        let w = c.sub(p);
1132        let (g, slope) = (w.dot(c1), c1.dot(c1) + w.dot(c2));
1133        if !(g.is_finite() && slope.is_finite()) || slope == 0.0 {
1134            wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: slope {slope:e} g {g:e} at t {t}"));
1135            return Solve::Unread;
1136        }
1137        if g.abs() <= 16.0 * f64::EPSILON * (p.length() + c.length()) * c1.length() {
1138            converged = true;
1139            break;
1140        }
1141        let step = g / slope;
1142        let before = t;
1143        t -= step;
1144        if !(t >= low && t <= high) {
1145            wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: window exit t {before} -> {t} (window [{low}, {high}])"));
1146            return if t.is_finite() { Solve::Window(t) } else { Solve::Unread };
1147        }
1148        // The uniform reader's stop, 1e-14 · (1 + |s|) in node units from
1149        // the stencil's base, read in LOCAL node units here (one node = half
1150        // the window), so it does not depend on where the native t's origin
1151        // is or on its scale.
1152        let node = 0.5 * (high - low);
1153        if (step / node).abs() <= 1e-14 * (1.0 + ((t - t0) / node).abs()) {
1154            converged = true;
1155            break;
1156        }
1157    }
1158    if !converged {
1159        wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: no convergence in 50 steps, t {t}"));
1160        return Solve::Unread;
1161    }
1162    let distance = curve(t).0.sub(p).length() - radius;
1163    if distance.is_finite() {
1164        Solve::Read(distance)
1165    } else {
1166        wall_decline(|| format!("native stencil {n} nodes base {base} j {j}: non-finite distance at t {t}"));
1167        Solve::Unread
1168    }
1169}
1170
1171/// The intervals the construction request refines (and, when it cannot,
1172/// reports): every interval whose DECLARED wall reading — as
1173/// `wall_interior_misses` returns it, `(left, declared, reading,
1174/// uncertainty)` — exceeds `model`. The one decision the request makes.
1175pub(in crate::blend) fn wall_model_missing(readings: &[(usize, f64, f64, f64)], model: f64) -> Vec<usize> {
1176    readings.iter().filter(|&&(_, declared, _, _)| !(declared <= model)).map(|&(left, _, _, _)| left).collect()
1177}
1178
1179/// Every fillet-wall interval's worst distance from the rolling ball's sweep,
1180/// read on the BUILT wall at the even fractions of the centre stencil
1181/// j/`subdivisions` (the centres at twice the probe density) and section
1182/// angles k/`probes` ([`declared_wall_probe`]), keyed by the interval's left
1183/// sample: each readable interval's largest DECLARED reading, largest reading
1184/// and largest uncertainty, and separately
1185/// every interval with an unreadable reading or a contact the stencil did not
1186/// solve. The wall shares the rails' `u` (one interpolation parameter for its
1187/// three rows).
1188#[allow(clippy::too_many_arguments)]
1189pub(super) fn wall_interior_misses(
1190    contacts: &crate::blend::chain::DenseContacts,
1191    samples: &[ChainSample],
1192    surface: &crate::NurbsSurface,
1193    low: f64,
1194    range: f64,
1195    radius: f64,
1196    subdivisions: usize,
1197    probes: usize,
1198    wrap: bool,
1199) -> (Vec<(usize, f64, f64, f64)>, Vec<usize>) {
1200    use std::collections::BTreeMap;
1201    // Each interval's centres at fractions 0, 1/k, …, 1: its stations at the
1202    // ends, the stencil's solved balls between.
1203    let mut nodes: BTreeMap<usize, (usize, Vec<Option<Vec3>>)> = BTreeMap::new();
1204    for (index, (left, next, fraction, _, _)) in contacts.solved.iter().enumerate() {
1205        let entry = nodes.entry(*left).or_insert_with(|| {
1206            let mut row = vec![None; subdivisions + 1];
1207            row[0] = Some(samples[*left].station.center);
1208            row[subdivisions] = Some(samples[*next].station.center);
1209            (*next, row)
1210        });
1211        let step = (fraction * subdivisions as f64).round() as usize;
1212        if (1..subdivisions).contains(&step) {
1213            entry.1[step] = contacts.centers.get(index).copied().flatten();
1214        }
1215    }
1216    let mut readable = Vec::new();
1217    let mut unreadable: Vec<usize> = contacts.unsolved_at.clone();
1218    for (&left, (next, row)) in &nodes {
1219        let Some(centres) = row.iter().copied().collect::<Option<Vec<Vec3>>>() else {
1220            unreadable.push(left);
1221            continue;
1222        };
1223        let (from, to) = (samples[left].parameter, samples[*next].parameter);
1224        let (mut declared, mut worst, mut spread) = (0.0_f64, 0.0_f64, 0.0_f64);
1225        let mut read = true;
1226        let prepared = PreparedWall::for_interval(&centres);
1227        'probes: for step in (2..subdivisions).step_by(2) {
1228            let global = from + (step as f64 / subdivisions as f64) * (to - from);
1229            let global = if wrap { global.rem_euclid(1.0) } else { global };
1230            let u = (global - low) / range;
1231            for k in 1..probes {
1232                let v = k as f64 / probes as f64;
1233                let reading = surface.evaluate(u, v).ok().and_then(|p| {
1234                    match &prepared {
1235                        Some(cache) => cache.read(p, step, step as f64 / subdivisions as f64, radius),
1236                        None => declared_wall_probe_located(p, &centres, step, step as f64 / subdivisions as f64, radius),
1237                    }
1238                });
1239                match reading {
1240                    Some((distance, gap)) => {
1241                        declared = declared.max(distance + gap);
1242                        worst = worst.max(distance);
1243                        spread = spread.max(gap);
1244                    }
1245                    None => {
1246                        read = false;
1247                        break 'probes;
1248                    }
1249                }
1250            }
1251        }
1252        if read {
1253            readable.push((left, declared, worst, spread));
1254        } else {
1255            unreadable.push(left);
1256        }
1257    }
1258    unreadable.sort_unstable();
1259    unreadable.dedup();
1260    (readable, unreadable)
1261}
1262
1263
1264fn build_closed_smooth_chain(
1265    solid: &BrepSolid,
1266    segments: &[ChainSegment<'_>],
1267    radius: f64,
1268    chamfer: bool,
1269    name: Option<&str>,
1270    budget: &ChainBudget,
1271    fold: FoldPolicy,
1272    bar: f64,
1273) -> Result<Rung, KernelRefusal> {
1274    let samples = match march_chain_typed(segments, radius, false, budget, fold) {
1275        Ok(samples) => samples,
1276        Err(ChainMarchError::Station(failure)) => return Ok(Rung::StationFailed(failure)),
1277        // Only local refinement reports incoherence; the march never does, so
1278        // this is its refusal exactly as `march_chain` would return it.
1279        Err(ChainMarchError::Incoherent(failure)) => return Err(failure.error),
1280        Err(ChainMarchError::Other(error)) => return Err(error),
1281    };
1282    CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(0));
1283    build_closed_chain_from_samples(solid, segments, radius, chamfer, name, budget, fold, bar, samples, 0)
1284}
1285
1286
1287thread_local! {
1288    /// The deepest local round any rung of the current ladder step attempted
1289    /// (reset as each rung's march is built from), so a rejected refinement
1290    /// reports the rounds it actually spent.
1291    static CHAIN_DEEPEST_ROUND: std::cell::Cell<usize> = const { std::cell::Cell::new(0) };
1292}
1293
1294// BREP lane L: begin lcount234_chain_local_continue_20261016
1295/// Whether a rung's local refinement may run past [`REFINE_ROUNDS`] while it
1296/// still improves (`BREP_CHAIN_LOCAL_CONTINUE=1`; off by default until its
1297/// census promotes it). The ladder repeats its local rounds from scratch on
1298/// every rung.
1299fn local_continue_on() -> bool {
1300    std::env::var("BREP_CHAIN_LOCAL_CONTINUE").as_deref() == Ok("1")
1301}
1302
1303/// The most local rounds a rung may take with the switch on.
1304const LOCAL_CONTINUE_MAX_ROUNDS: usize = 4 * REFINE_ROUNDS;
1305
1306/// Rounds without a new best stencil reading before the rung is spent.
1307const LOCAL_CONTINUE_STALL_ROUNDS: usize = 3;
1308
1309thread_local! {
1310    /// The best stencil reading of the rung being refined, and its round.
1311    static CHAIN_LOCAL_BEST: std::cell::Cell<(f64, usize)> = const { std::cell::Cell::new((f64::INFINITY, 0)) };
1312}
1313
1314/// Whether a rung's local rounds are spent with the switch on: past
1315/// [`LOCAL_CONTINUE_MAX_ROUNDS`], or [`LOCAL_CONTINUE_STALL_ROUNDS`] rounds
1316/// since the best stencil reading. Counted in rounds, never in time. Round 0
1317/// reads the halfway point only, a different instrument from the dense stencil
1318/// every later round reads, so the best is kept from round 1 on.
1319fn local_rounds_spent(round: usize, deviation: f64) -> bool {
1320    if round == 0 {
1321        CHAIN_LOCAL_BEST.with(|best| best.set((f64::INFINITY, 0)));
1322        return false;
1323    }
1324    let (_, best_round) = CHAIN_LOCAL_BEST.with(|best| {
1325        let (reading, at) = best.get();
1326        if deviation < reading {
1327            best.set((deviation, round));
1328            (deviation, round)
1329        } else {
1330            (reading, at)
1331        }
1332    });
1333    round >= LOCAL_CONTINUE_MAX_ROUNDS || round - best_round >= LOCAL_CONTINUE_STALL_ROUNDS
1334}
1335// BREP lane L: end lcount234_chain_local_continue_20261016
1336
1337/// Fit, measure and build one rung from its stations. A plain (non-carve)
1338/// rung whose fitted rails miss the halfway contacts by more than `bar`
1339/// first REFINES LOCALLY: one continuation station into each failing interval
1340/// (`march::insert_chain_stations`), up to [`REFINE_ROUNDS`] rounds and never
1341/// past [`CHAIN_CARVE_MAX_PER_SEGMENT`] stations in a segment, then measures
1342/// again. Only a rung whose local rounds are spent is too coarse, and the
1343/// ladder climbs as it always has. A rung that passes at once is built from
1344/// exactly the stations it was built from before.
1345#[allow(clippy::too_many_arguments)]
1346fn build_closed_chain_from_samples(
1347    solid: &BrepSolid,
1348    segments: &[ChainSegment<'_>],
1349    radius: f64,
1350    chamfer: bool,
1351    name: Option<&str>,
1352    budget: &ChainBudget,
1353    fold: FoldPolicy,
1354    bar: f64,
1355    mut samples: Vec<ChainSample>,
1356    round: usize,
1357) -> Result<Rung, KernelRefusal> {
1358    // The rung, for the traces and the typed reports; each segment's own
1359    // count is the budget's.
1360    let per_segment = budget.rung;
1361    CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round)));
1362
1363    // Global rows from the in-segment samples (wrapped-overlap closed fit).
1364    let mut global: Vec<(f64, &ChainSample)> = samples
1365        .iter()
1366        .filter(|sample| (0..budget.count(sample.segment) as isize).contains(&sample.position))
1367        .map(|sample| (sample.parameter.rem_euclid(1.0), sample))
1368        .collect();
1369    global.sort_by(|a, b| a.0.total_cmp(&b.0));
1370    let count = global.len();
1371    let degree = FIT_DEGREE;
1372    let mut extended_params = Vec::new();
1373    let mut samples_cr = Vec::new();
1374    let mut samples_mid = Vec::new();
1375    let mut samples_cs = Vec::new();
1376    let mut push = |station: &Station, parameter: f64| {
1377        extended_params.push(parameter);
1378        samples_cr.push(Vec4::from_point(station.p1, 1.0));
1379        samples_cs.push(Vec4::from_point(station.p2, 1.0));
1380        samples_mid.push(Vec4 {
1381            x: station.apex.x * station.weight,
1382            y: station.apex.y * station.weight,
1383            z: station.apex.z * station.weight,
1384            w: station.weight,
1385        });
1386    };
1387    for offset in (1..=degree).rev() {
1388        let (parameter, sample) = &global[count - offset];
1389        push(&sample.station, parameter - 1.0);
1390    }
1391    for (parameter, sample) in &global {
1392        push(&sample.station, *parameter);
1393    }
1394    // Close the loop: repeat the first station at parameter 1, then the
1395    // wrap continuation.
1396    push(&global[0].1.station, global[0].0 + 1.0);
1397    for offset in 1..=degree {
1398        let (parameter, sample) = &global[offset];
1399        push(&sample.station, parameter + 1.0);
1400    }
1401    let low = extended_params[0];
1402    let high = *extended_params.last().unwrap();
1403    let range = high - low;
1404    let normalized: Vec<f64> = extended_params
1405        .iter()
1406        .map(|parameter| (parameter - low) / range)
1407        .collect();
1408    let seam_low = (0.0 - low) / range;
1409    let seam_high = (1.0 - low) / range;
1410    let fit_row = |row: &[Vec4]| -> Result<NurbsCurve, KernelRefusal> {
1411        let curve = fit::interpolate_homogeneous(row, degree, &normalized).or_refuse(KernelStage::Refine, "interpolate_homogeneous")?;
1412        let (_, tail) = curve.split(seam_low).or_refuse(KernelStage::Refine, "split")?;
1413        let (middle, _) = tail.split(seam_high).or_refuse(KernelStage::Refine, "split")?;
1414        Ok(middle)
1415    };
1416    let cr = fit_row(&samples_cr)?;
1417    let cs = fit_row(&samples_cs)?;
1418    let mid = if chamfer {
1419        None
1420    } else {
1421        Some(fit_row(&samples_mid)?)
1422    };
1423    let u_domain = cr.domain().or_refuse(KernelStage::Refine, "domain")?;
1424    let surface = crate::blend::rows::surface_from_rows(degree, &cr, &cs, mid.as_ref(), true)?;
1425
1426    // The CARVE.  A wall that folds carries a lens of parameter inside the
1427    // swept ball volume; the blend is that wall with the lens cut away along
1428    // the crease where the two sheets meet, and the crease becomes an edge of
1429    // the blend with two coedges on this one face.  Traced on the surface that
1430    // was just fitted — not on the ideal envelope — so the two kept sheets
1431    // meet EXACTLY on the curve the trim is built from (`blend/carve.rs`).
1432    // EVERY rung is measured before anything is built on it, carve or not: the
1433    // rails are interpolants through the stations, and between stations they
1434    // can leave the rolling ball's own contacts — and so the carriers they are
1435    // supposed to ride — by far more than the kernel's intersection-fit
1436    // contract. A carve pays for surgery only on an accurate rung; a plain
1437    // collar pays nothing more than the measurement unless it misses.
1438    let (halfway, halfway_failing) = rail_deviation(&samples, &cr, &cs, low, range, budget, bar)?;
1439    // A rung that local insertion has refined is judged on the dense stencil of
1440    // exact contacts, not the halfway point alone (the NX collar's rung 24 read
1441    // 1.943844e-6 halfway and 2.009538e-6 at 7/16 against 2e-6). A rung that
1442    // passes without insertion is read exactly as it always was.
1443    // A FILLET wall is also judged between its rails, on every rung (carve
1444    // included): the rails can pass while the wall's interior — its weighted
1445    // apex row interpolated between stations — stands off the rolling ball
1446    // (the 20-degree seam-split row's accepted wall read 2.4335e-6 on a 16 x 9
1447    // probe grid and 2.4430e-6 at its worst interval's midpoint, mid-section,
1448    // against 2e-6). A chamfer is not a rolling ball and is not judged
1449    // against one. One stencil of exact contacts serves rails and wall.
1450    let stencil = if round > 0 {
1451        Some(crate::blend::chain::chain_dense_contacts(segments, radius, budget, &samples, LOCAL_RAIL_STENCIL)?)
1452    } else {
1453        None
1454    };
1455    // The wall's own centres, at twice the rail stencil's density: its probes
1456    // sit at the rail stencil's fractions, and each is read at both densities
1457    // (`declared_wall_probe`). The rails keep their own stencil, read exactly
1458    // as before.
1459    let wall_stencil = if !chamfer && wall_verdict_on() {
1460        Some(crate::blend::chain::chain_dense_contacts(segments, radius, budget, &samples, 2 * LOCAL_RAIL_STENCIL)?)
1461    } else {
1462        None
1463    };
1464    let (mut deviation, mut failing, _stencil_reads) = if round > 0 {
1465        let contacts = stencil.as_ref().expect("the local path solves the stencil");
1466        let (dense, unreadable) = dense_rail_misses(contacts, &samples, &cr, &cs, low, range, true)?;
1467        let reads = dense.len() + unreadable.len() + contacts.unsolved;
1468        let failed_to_read: Vec<usize> = contacts.unsolved_at.iter().chain(&unreadable).copied().collect();
1469        let verdict = combine_rail_verdict((halfway, halfway_failing), &dense, &failed_to_read, bar);
1470        crate::blend::carve::carve_trace(format_args!(
1471            "blend chain local rails: rung {per_segment} round {round}: stencil verdict {:.9e} (halfway {halfway:.9e}), \
1472             {} failing interval(s), {} unsolved contact(s)",
1473            verdict.0, verdict.1.len(), contacts.unsolved
1474        ));
1475        (verdict.0, verdict.1, reads)
1476    } else {
1477        (halfway, halfway_failing, 0)
1478    };
1479    // The wall between the rails, joined FAIL-CLOSED: an interval whose wall
1480    // reads over `bar`, or cannot be read, fails exactly as a rail interval
1481    // does, and the ladder (local insertion, climb, the existing refusal)
1482    // decides as it always has. Nothing else about the rung changes: a wall
1483    // inside `bar` adds reads, never stations.
1484    let mut _wall_worst = 0.0_f64;
1485    let mut _wall_spread = 0.0_f64;
1486    let mut wall_readings: Vec<(usize, f64, f64, f64)> = Vec::new();
1487    if let Some(contacts) = wall_stencil.as_ref() {
1488        let (wall, unreadable) =
1489            wall_interior_misses(contacts, &samples, &surface, low, range, radius, 2 * LOCAL_RAIL_STENCIL, WALL_SECTION_PROBES, true);
1490        wall_readings = wall.clone();
1491        // Acceptance on the DECLARED reading against the unchanged bar.
1492        for &(left, declared, _, uncertainty) in &wall {
1493            if declared > bar {
1494                failing.push(left);
1495            }
1496            _wall_worst = _wall_worst.max(declared);
1497            _wall_spread = _wall_spread.max(uncertainty);
1498        }
1499        failing.extend_from_slice(&unreadable);
1500        failing.sort_unstable();
1501        failing.dedup();
1502        deviation = deviation.max(_wall_worst);
1503        crate::blend::carve::carve_trace(format_args!(
1504            "blend chain wall: rung {per_segment} round {round}: wall interior worst {_wall_worst:.9e} over {} \
1505             interval(s), {} unreadable, bar {bar:.3e}",
1506            wall.len(),
1507            unreadable.len()
1508        ));
1509    }
1510    if deviation > bar || !failing.is_empty() {
1511        let in_segment = |segment: usize| {
1512            samples
1513                .iter()
1514                .filter(|sample| sample.segment == segment && (0..budget.count(segment) as isize).contains(&sample.position))
1515                .count()
1516        };
1517        let within_ceiling = (0..segments.len()).all(|segment| {
1518            in_segment(segment) + failing.iter().filter(|&&index| samples[index].segment == segment).count()
1519                <= CHAIN_CARVE_MAX_PER_SEGMENT
1520        });
1521        crate::blend::carve::carve_trace(format_args!(
1522            "blend chain local rails: rung {per_segment} round {round} stations {} failing {} \
1523             worst {deviation:.9e} bar {bar:.9e} ceiling {}",
1524            samples.len(),
1525            failing.len(),
1526            if within_ceiling { "ok" } else { "reached" }
1527        ));
1528        // BREP lane L: begin lcount234_chain_local_continue_20261016
1529        let rounds_spent = if local_continue_on() {
1530            local_rounds_spent(round, deviation)
1531        } else {
1532            round >= REFINE_ROUNDS
1533        };
1534        // BREP lane L: end lcount234_chain_local_continue_20261016
1535        if fold != FoldPolicy::Refuse || rounds_spent || !within_ceiling {
1536            return Ok(Rung::TooCoarse { deviation });
1537        }
1538        // The insertion is the next round's attempt, whatever it returns.
1539        CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round + 1)));
1540        match insert_chain_stations(segments, radius, budget, &mut samples, &failing) {
1541            Ok(()) => {}
1542            Err(ChainMarchError::Station(failure)) => return Ok(Rung::StationFailed(failure)),
1543            // An inserted station that does not reach its retained right
1544            // neighbour's branch: this rung cannot be refined coherently, so it
1545            // is as coarse as its measured rails say, and the ladder decides.
1546            Err(ChainMarchError::Incoherent(failure)) => {
1547                crate::blend::carve::carve_trace(format_args!(
1548                    "blend chain local rails: rung {per_segment} round {round}: segment {} station \
1549                     at t {:.9} does not reach its retained right neighbour: {}",
1550                    failure.segment, failure.t, failure.error
1551                ));
1552                return Ok(Rung::TooCoarse { deviation });
1553            }
1554            Err(ChainMarchError::Other(error)) => return Err(error),
1555        }
1556        assign_chain_parameters(&mut samples, segments.len(), budget, false);
1557        return build_closed_chain_from_samples(
1558            solid, segments, radius, chamfer, name, budget, fold, bar, samples, round + 1,
1559        );
1560    }
1561    crate::blend::carve::carve_trace(format_args!(
1562        "blend chain: {per_segment} stations per segment leaves the rails {deviation:.6e} \
1563         from the rolling ball, inside {bar:.3e}"
1564    ));
1565    // The CONSTRUCTION request, plain fillet rungs only: an ACCEPTED wall is
1566    // refined toward `KernelTolerances::model` — the identity tolerance, the
1567    // same 1e-7 the pcurve construction floor refines to under its looser
1568    // acceptance — locally, in exactly the intervals whose DECLARED reading
1569    // (`declared_wall_probe`: the proxy reading plus its uncertainty from
1570    // stencil spread AND change with node density) is over it: one rule for
1571    // the acceptance verdict, this request and its report. An estimate with
1572    // its uncertainty, NOT a proven bound.
1573    // Acceptance is unchanged: the refined rung must pass it again, and every
1574    // retained station, branch guard and ceiling holds. A wall already inside
1575    // `model` everywhere inserts nothing and builds exactly as it did. When
1576    // refinement cannot run or does not build, THIS accepted rung stands and
1577    // says so, typed (`WallModelReport`); never silently.
1578    let model = crate::KernelTolerances::for_solid(solid, 1e-7).model;
1579    let stations_now = (0..segments.len())
1580        .map(|segment| {
1581            samples
1582                .iter()
1583                .filter(|sample| sample.segment == segment && (0..budget.count(segment) as isize).contains(&sample.position))
1584                .count()
1585        })
1586        .max()
1587        .unwrap_or(0);
1588    let mut unmet: Option<(f64, crate::BudgetReason)> = None;
1589    // What the request found short, named in its report: the detail, which
1590    // component the residual is the reading of, and how many intervals it
1591    // could not read (then the true worst is unknown).
1592    let mut request_detail: Option<String> = None;
1593    let mut request_component = crate::MeasuredComponent::Wall;
1594    let mut request_unread = 0usize;
1595    if !chamfer && fold == FoldPolicy::Refuse && model_request_on() {
1596        let wall_missing = wall_model_missing(&wall_readings, model);
1597        // The RAILS against the same request: the shipped support curves ARE
1598        // these rails (`chain_surgery` splits `cr`/`cs` into the support
1599        // edges), so a rail standing off the exact construction ships that
1600        // miss on its carrier. Read against the exact contacts the wall
1601        // reading already solved (j/32 of every interval); an interval whose
1602        // rail reads over `model`, or cannot be read, is short exactly as a
1603        // wall interval is. Acceptance (`bar`) is unchanged.
1604        let (rail_dense, rail_unread) = match wall_stencil.as_ref() {
1605            Some(contacts) => dense_rail_misses(contacts, &samples, &cr, &cs, low, range, true)?,
1606            None => (Vec::new(), Vec::new()),
1607        };
1608        // Per contact, keyed by the interval's left sample: distinct intervals.
1609        let rail_over: std::collections::BTreeSet<usize> =
1610            rail_dense.iter().filter(|&&(_, miss)| !(miss <= model)).map(|&(left, _)| left).collect();
1611        let rail_unread_intervals: std::collections::BTreeSet<usize> = rail_unread.iter().copied().collect();
1612        let rail_missing: Vec<usize> = rail_over.union(&rail_unread_intervals).copied().collect();
1613        let mut missing: Vec<usize> = wall_missing.iter().chain(&rail_missing).copied().collect();
1614        missing.sort_unstable();
1615        missing.dedup();
1616        if !missing.is_empty() {
1617            // The same declared reading, over the intervals short.
1618            let wall_short = wall_readings
1619                .iter()
1620                .filter(|reading| wall_missing.contains(&reading.0))
1621                .map(|&(_, declared, _, _)| declared)
1622                .fold(0.0_f64, f64::max);
1623            let rail_short = rail_dense.iter().filter(|&&(_, miss)| !(miss <= model)).map(|&(_, miss)| miss).fold(0.0_f64, f64::max);
1624            let short = wall_short.max(rail_short);
1625            request_component = if rail_short > wall_short { crate::MeasuredComponent::Rails } else { crate::MeasuredComponent::Wall };
1626            request_unread = rail_unread_intervals.len();
1627            request_detail = Some(format!(
1628                "short of {model:.3e}: rails {rail_short:.3e} off the exact contacts over {} interval(s) ({} unread), \
1629                 wall {wall_short:.3e} over {} interval(s)",
1630                rail_over.len(),
1631                rail_unread_intervals.len(),
1632                wall_missing.len()
1633            ));
1634            let within_ceiling = (0..segments.len()).all(|segment| {
1635                let present = samples
1636                    .iter()
1637                    .filter(|sample| sample.segment == segment && (0..budget.count(segment) as isize).contains(&sample.position))
1638                    .count();
1639                present + missing.iter().filter(|&&index| samples[index].segment == segment).count()
1640                    <= CHAIN_CARVE_MAX_PER_SEGMENT
1641            });
1642            crate::blend::carve::carve_trace(format_args!(
1643                "blend chain wall model: rung {per_segment} round {round}: {} interval(s) short of {model:.3e} \
1644                 (worst reading + spread {short:.3e}), ceiling {}",
1645                missing.len(),
1646                if within_ceiling { "ok" } else { "reached" }
1647            ));
1648            unmet = if round >= model_rounds_limit() {
1649                Some((short, crate::BudgetReason::RoundsSpent))
1650            } else if !within_ceiling {
1651                Some((short, crate::BudgetReason::StationCeiling))
1652            } else {
1653                let mut refined = samples.clone();
1654                // The insertion is the next round's attempt, whatever it returns.
1655                CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round + 1)));
1656                match insert_chain_stations(segments, radius, budget, &mut refined, &missing) {
1657                    Ok(()) => {
1658                        assign_chain_parameters(&mut refined, segments.len(), budget, false);
1659                        match build_closed_chain_from_samples(
1660                            solid, segments, radius, chamfer, name, budget, fold, bar, refined, round + 1,
1661                        ) {
1662                            Ok(Rung::Built(built)) => return Ok(Rung::Built(built)),
1663                            _ => Some((short, crate::BudgetReason::Rejected)),
1664                        }
1665                    }
1666                    Err(_) => {
1667                        // The insertion itself was the attempt.
1668                        CHAIN_DEEPEST_ROUND.with(|deepest| deepest.set(deepest.get().max(round + 1)));
1669                        Some((short, crate::BudgetReason::Incoherent))
1670                    }
1671                }
1672            };
1673        }
1674    }
1675    // `BREP_DEBUG_CHAIN_DENSE_RAILS=k` (diagnostic only): the accepted rung's
1676    // rails against exact contacts at k-1 interior fractions of every interval,
1677    // where `rail_deviation` reads only the halfway one. Nothing reads it to
1678    // decide.
1679    if let Some(subdivisions) = std::env::var("BREP_DEBUG_CHAIN_DENSE_RAILS").ok().and_then(|value| value.parse::<usize>().ok()) {
1680        let contacts = crate::blend::chain::chain_dense_contacts(segments, radius, budget, &samples, subdivisions.max(2))?;
1681        let mut worst = (0.0_f64, 0usize, 0isize, 0.0_f64);
1682        let mut unreadable = 0usize;
1683        for (left, next, fraction, p1, p2) in &contacts.solved {
1684            let (from, to) = (samples[*left].parameter, samples[*next].parameter);
1685            let global = (from + fraction * (to - from)).rem_euclid(1.0);
1686            let seed = (global - low) / range;
1687            let reach = (to - from).abs() / range;
1688            match rail_pair_miss(closest_on_rail(&cr, *p1, seed, reach).unwrap_or(f64::NAN), closest_on_rail(&cs, *p2, seed, reach).unwrap_or(f64::NAN)) {
1689                Some(miss) if miss > worst.0 => worst = (miss, samples[*left].segment, samples[*left].position, *fraction),
1690                Some(_) => {}
1691                None => unreadable += 1,
1692            }
1693        }
1694        eprintln!(
1695            "blend chain dense rails: rung {per_segment} round {round}: {} contacts at {subdivisions} per interval, \
1696             worst {:.9e} (halfway reading {halfway:.9e}, bar {bar:.3e}) at segment {} position {} fraction {:.4}; {} unsolved, {unreadable} unreadable",
1697            contacts.solved.len(), worst.0, worst.1, worst.2, worst.3, contacts.unsolved
1698        );
1699    }
1700    let crease = if fold == FoldPolicy::Carve {
1701        // THE BAR IS THE KERNEL'S OWN CONTRACT FOR THE QUANTITY MEASURED. What
1702        // the fit is held to is "does this coedge's curve-on-surface, pushed
1703        // back to 3-D, still trace the same locus as the edge's 3-D curve?",
1704        // which is `pcurve_consistency` — not `intersection_fit`, which is the
1705        // accuracy of a fitted intersection CURVE and is what the re-marched
1706        // wall itself is held to (`carve_ladder`). Measured on the
1707        // 2026-09-02 collar the crease's pcurves reach 5.4e-6 against that
1708        // 4e-3, so the margin is three decades and the bar is not what decides
1709        // the case either way.
1710        let bar = crate::KernelTolerances::for_solid(solid, 1e-7).pcurve_consistency;
1711        let Some(carved) = crate::blend::carve::trace_wall_crease(&surface, radius.abs())? else {
1712            return Err(KernelRefusal::internal(KernelStage::Refine, "carve_no_fold", format!(
1713                "blend carve: the wall re-marched at {per_segment} stations per segment carries                  no fold to carve, but the march at {CHAIN_PER_SEGMENT} refused one"
1714            )));
1715        };
1716        Some(crate::blend::carve::fit_crease(&surface, &carved, bar)?)
1717    } else {
1718        None
1719    };
1720
1721    // Single-face sides (e.g. one cap around the whole rim) get ONE uv
1722    // fit in the rows' own normalized space — identical parameterization
1723    // to cr/cs, so the closed support piece and its pcurve agree exactly.
1724    let single_face = |side: usize| -> bool {
1725        let first_id = segments[0].mate(side).face.id;
1726        segments
1727            .iter()
1728            .all(|segment| segment.mate(side).face.id == first_id)
1729    };
1730    let mut whole_pcurves: [Option<NurbsCurve>; 2] = [None, None];
1731    for side in 0..2 {
1732        if !single_face(side) {
1733            continue;
1734        }
1735        let mut row = Vec::new();
1736        let select = |station: &Station| -> [f64; 2] {
1737            if side == 0 {
1738                station.uv1
1739            } else {
1740                station.uv2
1741            }
1742        };
1743        let mut push_uv = |station: &Station| {
1744            let uv = select(station);
1745            row.push(Vec4::from_point(Vec3::new(uv[0], uv[1], 0.0), 1.0));
1746        };
1747        for offset in (1..=degree).rev() {
1748            push_uv(&global[count - offset].1.station);
1749        }
1750        for (_, sample) in &global {
1751            push_uv(&sample.station);
1752        }
1753        push_uv(&global[0].1.station);
1754        for offset in 1..=degree {
1755            push_uv(&global[offset].1.station);
1756        }
1757        whole_pcurves[side] = Some(fit_row(&row)?);
1758    }
1759
1760    // Per-face pcurve fits over each segment's full sample window
1761    // (in-segment + overshoot), in global parameters.
1762    let mut pcurves1 = Vec::with_capacity(segments.len());
1763    let mut pcurves2 = Vec::with_capacity(segments.len());
1764    for segment in 0..segments.len() {
1765        let mut window: Vec<&ChainSample> = samples
1766            .iter()
1767            .filter(|sample| sample.segment == segment)
1768            .collect();
1769        window.sort_by(|a, b| a.parameter.total_cmp(&b.parameter));
1770        let params: Vec<f64> = window.iter().map(|sample| sample.parameter).collect();
1771        let low = params[0];
1772        let high = *params.last().unwrap();
1773        let normalized: Vec<f64> = params
1774            .iter()
1775            .map(|parameter| (parameter - low) / (high - low))
1776            .collect();
1777        let fit_uv = |select: &dyn Fn(&Station) -> [f64; 2]| -> Result<NurbsCurve, KernelRefusal> {
1778            let points: Vec<Vec4> = window
1779                .iter()
1780                .map(|sample| {
1781                    let uv = select(&sample.station);
1782                    Vec4::from_point(Vec3::new(uv[0], uv[1], 0.0), 1.0)
1783                })
1784                .collect();
1785            let curve = fit::interpolate_homogeneous(&points, degree, &normalized).or_refuse(KernelStage::Refine, "interpolate_homogeneous")?;
1786            // Re-express in GLOBAL parameters: the curve's [0,1] domain
1787            // corresponds to [low, high]; keep as-is and remember the
1788            // affine map through the window bounds.
1789            Ok(curve)
1790        };
1791        pcurves1.push((low, high, fit_uv(&|station| station.uv1)?));
1792        pcurves2.push((low, high, fit_uv(&|station| station.uv2)?));
1793    }
1794
1795    let wall_surface = unmet.as_ref().map(|_| surface.clone());
1796    let built = chain_surgery(
1797        solid,
1798        segments,
1799        ChainRows {
1800            surface,
1801            cr,
1802            cs,
1803            u_domain,
1804            fit_low: low,
1805            fit_range: range,
1806            pcurves1,
1807            pcurves2,
1808            whole_pcurves,
1809            crease,
1810        },
1811        name,
1812    )?;
1813    if let (Some((residual, reason)), Some(surface)) = (unmet, wall_surface) {
1814        crate::blend::record_wall_model_report(crate::blend::WallModelReport {
1815            name: name.map(str::to_string),
1816            surface,
1817            request: model,
1818            residual,
1819            detail: request_detail,
1820            budget: crate::ApproximationBudget {
1821                reason,
1822                // Rounds ATTEMPTED (the deepest any refinement of this rung
1823                // reached); the stations and residual are this shipped rung's.
1824                rounds_used: CHAIN_DEEPEST_ROUND.with(|deepest| deepest.get()).max(round),
1825                rounds_limit: model_rounds_limit(),
1826                stations: stations_now,
1827                station_limit: CHAIN_CARVE_MAX_PER_SEGMENT,
1828                mechanism: crate::BudgetMechanism::LocalRounds,
1829                measured_component: request_component,
1830                unread: request_unread,
1831            },
1832        });
1833    }
1834    Ok(Rung::Built(built))
1835}
1836
1837pub(super) struct ChainRows {
1838    pub(super) surface: NurbsSurface,
1839    pub(super) cr: NurbsCurve,
1840    pub(super) cs: NurbsCurve,
1841    pub(super) u_domain: [f64; 2],
1842    /// Affine map from the rows' fit space to GLOBAL chain parameters:
1843    /// global = fit_low + p * fit_range (the global fit normalised its
1844    /// wrap-extended parameters onto [0, 1] before interpolating).
1845    pub(super) fit_low: f64,
1846    pub(super) fit_range: f64,
1847    /// Per-segment (window_low, window_high, uv curve over [0,1]) in
1848    /// GLOBAL parameters.
1849    pub(super) pcurves1: Vec<(f64, f64, NurbsCurve)>,
1850    pub(super) pcurves2: Vec<(f64, f64, NurbsCurve)>,
1851    /// Whole-chain uv fits (rows' fit space) for single-face sides.
1852    pub(super) whole_pcurves: [Option<NurbsCurve>; 2],
1853    /// The crease a CARVED wall is cut along: the inner loop of the blend
1854    /// face, one edge carrying two coedges of that same face.
1855    pub(super) crease: Option<crate::blend::carve::CarvedCrease>,
1856}
1857
1858impl ChainRows {
1859    pub(super) fn to_global(&self, fit_parameter: f64) -> f64 {
1860        self.fit_low + fit_parameter * self.fit_range
1861    }
1862}
1863
1864/// Extract the pcurve portion for a global window [a, b] from a
1865/// segment's fitted uv curve, re-parameterised so its domain maps
1866/// affinely onto the portion (matching the split support edge).
1867pub(super) fn pcurve_portion(fitted: &(f64, f64, NurbsCurve), a: f64, b: f64) -> Result<NurbsCurve, KernelRefusal> {
1868    let (low, high, curve) = fitted;
1869    // The chain parameterisation wraps with period 1; pick the period
1870    // copy of the window that overlaps this segment's fit range.
1871    let to_local = |value: f64| {
1872        let mut best = value;
1873        for candidate in [value, value + 1.0, value - 1.0] {
1874            if candidate >= low - 1e-9 && candidate <= high + 1e-9 {
1875                best = candidate;
1876                break;
1877            }
1878        }
1879        ((best - low) / (high - low)).clamp(0.0, 1.0)
1880    };
1881    let mut local_a = to_local(a);
1882    let mut local_b = to_local(b);
1883    if local_a > local_b {
1884        std::mem::swap(&mut local_a, &mut local_b);
1885    }
1886    let epsilon = 1e-9;
1887    let (_, tail) = if local_a > epsilon {
1888        curve.split(local_a).or_refuse(KernelStage::Refine, "split")?
1889    } else {
1890        (curve.clone(), curve.clone())
1891    };
1892    let tail = if local_a > epsilon {
1893        tail
1894    } else {
1895        curve.clone()
1896    };
1897    let portion = if local_b < 1.0 - epsilon {
1898        tail.split(local_b).or_refuse(KernelStage::Refine, "split")?.0
1899    } else {
1900        tail
1901    };
1902    Ok(portion)
1903}
1904
1905/// The single boundary edge of a side's mate face(s) incident to a chain
1906/// vertex, EXCLUDING the two chain edges meeting there — i.e. the SEAM
1907/// (same face on both segments) or the SPOKE (face changes) that the
1908/// support curve must cross.  `None` when the chain simply flows through
1909/// (same face, no seam — e.g. a stadium cap rim vertex).
1910pub(super) fn cross_edge_at(
1911    solid: &BrepSolid,
1912    face_before: &FaceRecord,
1913    loop_before: usize,
1914    face_after: &FaceRecord,
1915    loop_after: usize,
1916    vertex: u64,
1917    before_edge: u64,
1918    after_edge: u64,
1919) -> Result<Option<u64>, KernelRefusal> {
1920    let mut found: Option<u64> = None;
1921    let mut consider = |face: &FaceRecord, loop_index: usize| -> Result<(), KernelRefusal> {
1922        for coedge in &face.loops[loop_index].coedges {
1923            if coedge.edge_id == before_edge || coedge.edge_id == after_edge {
1924                continue;
1925            }
1926            let edge = solid
1927                .edges
1928                .iter()
1929                .find(|edge| edge.id == coedge.edge_id)
1930                .ok_or(KernelRefusal::internal(KernelStage::Refine, "missing_edge", "blend: loop references missing edge"))?;
1931            if edge.start_vertex_id == vertex || edge.end_vertex_id == vertex {
1932                if found.is_some() && found != Some(edge.id) {
1933                    return Err(KernelRefusal::unsupported(KernelStage::Classify, "chain_cross_edges", "blend: multiple cross edges at a chain vertex"));
1934                }
1935                found = Some(edge.id);
1936            }
1937        }
1938        Ok(())
1939    };
1940    consider(face_before, loop_before)?;
1941    if !(face_after.id == face_before.id && loop_after == loop_before) {
1942        consider(face_after, loop_after)?;
1943    }
1944    Ok(found)
1945}
1946
1947/// The nearest periodic image of a seam boundary to `value` — the boundary a
1948/// crossing endpoint must be pinned to, read off its ADJACENT INTERIOR sample.
1949fn nearest_seam_image(value: f64, low: f64, high: f64) -> f64 {
1950    let period = high - low;
1951    low + ((value - low) / period).round() * period
1952}
1953
1954/// The mate-face pcurve of a support PIECE, built by projecting the piece's
1955/// 3D curve onto the (analytic) mate surface and fitting.  Because a chain
1956/// support rides EXACTLY on its mate carrier, the projection is exact; the
1957/// periodic track is unwrapped so it never jumps a period across the carrier
1958/// seam, and crossing endpoints are pinned ONTO the seam so the piece and the
1959/// trimmed seam edge close the loop in parameter space.
1960///
1961/// The carrier is closed in EITHER direction, and the crossed seam is not
1962/// always the u meridian.  A cylinder/cone/sphere is closed in u only, so its
1963/// seam is the u0 ≡ u1 meridian — but a TORUS is biperiodic, and a tube welded
1964/// through a plane crosses the torus' v0 ≡ v1 PARALLEL (its outer equator)
1965/// instead.  Pinning u there drags the endpoint a half-carrier away from its
1966/// own 3D vertex (the collar loop then closes on a uv corner that is nowhere
1967/// near the vertex, and the mate face's trim is garbage), so the pinned AXIS is
1968/// chosen GEOMETRICALLY: whichever closed direction's seam actually passes
1969/// through the endpoint.  A carrier closed in u only can only ever take the u
1970/// branch, exactly as before.
1971pub(in crate::blend) fn project_piece_pcurve(
1972    piece: &NurbsCurve,
1973    surface: &NurbsSurface,
1974    start_is_crossing: bool,
1975    end_is_crossing: bool,
1976) -> Result<NurbsCurve, KernelRefusal> {
1977    project_piece_pcurve_from(piece, surface, start_is_crossing, end_is_crossing, None)
1978}
1979
1980
1981/// A trim pcurve `old` as a GUIDE for [`project_piece_pcurve_from`]: at track
1982/// fraction `f` its uv at the same fraction of its own domain, closed
1983/// directions wrapped and open ones clamped into `surface`'s domain; `None`
1984/// where that uv is not finite.
1985pub(in crate::blend) fn pcurve_guide<'a>(old: &'a NurbsCurve, surface: &NurbsSurface) -> Result<impl Fn(f64) -> Option<[f64; 2]> + 'a, KernelRefusal> {
1986    let [o0, o1] = old.domain().or_refuse(KernelStage::Refine, "domain")?;
1987    let (closed_u, closed_v) = surface.closed_directions().or_refuse(KernelStage::Refine, "closed_directions")?;
1988    let [u0, u1] = surface.domain_u().or_refuse(KernelStage::Refine, "domain_u")?;
1989    let [v0, v1] = surface.domain_v().or_refuse(KernelStage::Refine, "domain_v")?;
1990    let into = move |value: f64, closed: bool, low: f64, high: f64| {
1991        if closed && high > low {
1992            low + (value - low).rem_euclid(high - low)
1993        } else {
1994            value.clamp(low, high)
1995        }
1996    };
1997    Ok(move |fraction: f64| {
1998        let uv = old.evaluate(o0 + (o1 - o0) * fraction).ok()?;
1999        let seed = [into(uv.x, closed_u, u0, u1), into(uv.y, closed_v, v0, v1)];
2000        seed.iter().all(|value| value.is_finite()).then_some(seed)
2001    })
2002}
2003
2004/// [`project_piece_pcurve`], every sample's foot sought from `guide` at its
2005/// fraction (a uv on the branch the piece is known to run on — a support
2006/// trim's own march pcurve) instead of from its neighbour's foot. A foot
2007/// continued from its neighbour stays wherever that neighbour's did: at a
2008/// STATIONARY edge of the carrier (an extrusion of a profile whose derivative
2009/// vanishes at its start, S_u = 0 along u = 0) a foot seeded on that edge has
2010/// no gradient to leave it by, and every later sample stayed there with it
2011/// (the stationary-start prism's refit: one foot, [0, 0.9], for the whole
2012/// rail, 6.0 off it at the far end). Without a guide, as before.
2013pub(in crate::blend) fn project_piece_pcurve_from(
2014    piece: &NurbsCurve,
2015    surface: &NurbsSurface,
2016    start_is_crossing: bool,
2017    end_is_crossing: bool,
2018    guide: Option<&dyn Fn(f64) -> Option<[f64; 2]>>,
2019) -> Result<NurbsCurve, KernelRefusal> {
2020    project_piece_pcurve_within(piece, surface, start_is_crossing, end_is_crossing, guide, &[])
2021}
2022
2023/// [`project_piece_pcurve_from`], its track first refined (two bisections)
2024/// in the dense interval holding each fraction of `refine_near` — where a
2025/// caller that read the SHIPPED curve again found it over the floor between
2026/// the fitter's own stencil points — and then fitted at the same floor.
2027/// The code of a support pcurve fit that misses the refinement floor.
2028pub(crate) const CHAIN_PCURVE_FLOOR_WHAT: &str = "chain_pcurve_floor";
2029
2030pub(in crate::blend) fn project_piece_pcurve_within(
2031    piece: &NurbsCurve,
2032    surface: &NurbsSurface,
2033    start_is_crossing: bool,
2034    end_is_crossing: bool,
2035    guide: Option<&dyn Fn(f64) -> Option<[f64; 2]>>,
2036    refine_near: &[f64],
2037) -> Result<NurbsCurve, KernelRefusal> {
2038    project_piece_pcurve_track(piece, surface, start_is_crossing, end_is_crossing, guide, refine_near, false)
2039}
2040
2041/// [`project_piece_pcurve_within`], retried where it misses the floor on a
2042/// track that also samples the piece's own knots (lane L, 2026-10-16). A track
2043/// sampled uniformly in the piece's parameter cannot follow a rail whose own
2044/// knot spans are finer than its spacing — a closed march crowds its stations
2045/// where the rim's curvature jumps (G1 not G2 at its seam: rail spans of
2046/// 3.6e-5 against a track spacing of 2.4e-4, the fit stalled 1.7e-5 off). Only
2047/// for a caller whose floor miss is TERMINAL: a caller that answers the miss
2048/// itself (a ladder that refines its march) keeps its own answer. A fit
2049/// already on the floor never takes the second track; a retry that misses
2050/// returns the uniform refusal unchanged.
2051pub(in crate::blend) fn project_piece_pcurve_within_knot_retry(
2052    piece: &NurbsCurve,
2053    surface: &NurbsSurface,
2054    start_is_crossing: bool,
2055    end_is_crossing: bool,
2056    guide: Option<&dyn Fn(f64) -> Option<[f64; 2]>>,
2057    refine_near: &[f64],
2058) -> Result<NurbsCurve, KernelRefusal> {
2059    match project_piece_pcurve_track(piece, surface, start_is_crossing, end_is_crossing, guide, refine_near, false) {
2060        Err(uniform) if matches!(&uniform.class, crate::RefusalClass::NonConvergence { what } if what == CHAIN_PCURVE_FLOOR_WHAT) => {
2061            project_piece_pcurve_track(piece, surface, start_is_crossing, end_is_crossing, guide, refine_near, true).map_err(|_| uniform)
2062        }
2063        other => other,
2064    }
2065}
2066
2067/// [`project_piece_pcurve_within`]'s track and fit, the track sampled
2068/// uniformly in the piece's parameter, and with `knot_aware` at the piece's
2069/// interior knots too.
2070#[allow(clippy::too_many_arguments)]
2071fn project_piece_pcurve_track(
2072    piece: &NurbsCurve,
2073    surface: &NurbsSurface,
2074    start_is_crossing: bool,
2075    end_is_crossing: bool,
2076    guide: Option<&dyn Fn(f64) -> Option<[f64; 2]>>,
2077    refine_near: &[f64],
2078    knot_aware: bool,
2079) -> Result<NurbsCurve, KernelRefusal> {
2080    let [u0, u1] = surface.domain_u().or_refuse(KernelStage::Refine, "domain_u")?;
2081    let [v0, v1] = surface.domain_v().or_refuse(KernelStage::Refine, "domain_v")?;
2082    let (_, closed_v) = surface.closed_directions().or_refuse(KernelStage::Refine, "closed_directions")?;
2083    let period = u1 - u0;
2084    let period_v = v1 - v0;
2085    let [d0, d1] = piece.domain().or_refuse(KernelStage::Refine, "domain")?;
2086    // Dense reference projection; the pcurve is fit from an adaptively
2087    // thinned subset kept as coarse as tolerance allows.
2088    //
2089    // The track is as dense as the PIECE is, not a fixed 96. A carved wall's
2090    // re-marched rail carries hundreds of control points (574 on the
2091    // 2026-09-02 two-torus collar), and a pcurve fitted and checked at 97
2092    // samples of it was accepted 1.19e-3 off the rail it claims.
2093    // The knot-aware retry's track holds the piece's knots on top of the
2094    // uniform samples; its repair may use twice the uniform ceiling.
2095    let ceiling = if knot_aware { 2 * PIECE_TRACK_CEILING } else { PIECE_TRACK_CEILING };
2096    let uniform = (8 * piece.control_points.len()).clamp(96, PIECE_TRACK_CEILING);
2097    let mut fractions: Vec<f64> = (0..=uniform).map(|section| section as f64 / uniform as f64).collect();
2098    if knot_aware {
2099        // The piece's distinct interior knots, as fractions, merged in where no
2100        // uniform sample already sits within a millionth of the span.
2101        let span = d1 - d0;
2102        let mut knots: Vec<f64> = piece.knots.iter().copied().filter(|knot| *knot > d0 && *knot < d1).map(|knot| (knot - d0) / span).collect();
2103        knots.dedup();
2104        fractions.extend(knots);
2105        fractions.sort_by(|a, b| a.total_cmp(b));
2106        fractions.dedup_by(|b, a| (*b - *a).abs() <= 1e-6 / uniform as f64);
2107    }
2108    let dense = fractions.len() - 1;
2109    crate::blend::record_count("count.track.piece_samples", dense as u64 + 1);
2110    let mut proj_u = Vec::with_capacity(dense + 1);
2111    let mut proj_v = Vec::with_capacity(dense + 1);
2112    let mut point3 = Vec::with_capacity(dense + 1);
2113    let mut params = Vec::with_capacity(dense + 1);
2114    let mut previous_u: Option<f64> = None;
2115    // The foot each sample continues from: the track is ONE branch of the
2116    // carrier, so after the first sample each inversion starts from its
2117    // neighbour's foot rather than asking for the nearest point of the whole
2118    // surface, which on a carrier that comes back near itself is another sheet.
2119    let mut previous_foot: Option<[f64; 2]> = None;
2120    let mut feet: Vec<[f64; 2]> = Vec::with_capacity(dense + 1);
2121    let project_timer = crate::blend::BlendTimer::start("timing.track.project_ms");
2122    for section in 0..=dense {
2123        let fraction = fractions[section];
2124        // The uniform track keeps its own expression, bit for bit (the
2125        // product before the quotient); only knot-aware samples use the
2126        // fraction directly.
2127        let t = if knot_aware { d0 + (d1 - d0) * fraction } else { d0 + (d1 - d0) * section as f64 / dense as f64 };
2128        let point = piece.evaluate(t).or_refuse(KernelStage::Refine, "evaluate")?;
2129        let guide = guide.filter(|_| {
2130            true
2131        });
2132        // The guide's seed and the neighbour's foot are both tried and the
2133        // NEARER foot kept: the seeded projector is local and returns a
2134        // stalled Newton's iterate without saying so, so neither seed is
2135        // trusted alone (a foot on its carrier is nearer than any stall).
2136        let seeds: Vec<[f64; 2]> = guide.and_then(|guide| guide(fraction)).into_iter().chain(previous_foot).collect();
2137        let mut projection: Option<crate::projection::SurfaceProjection> = None;
2138        for seed in &seeds {
2139            let found = crate::projection::project_point_to_surface_from_seed(surface, point, *seed).or_refuse(KernelStage::Refine, "project_point_to_surface_from_seed")?;
2140            if projection.as_ref().map_or(true, |known| found.distance < known.distance) {
2141                projection = Some(found);
2142            }
2143        }
2144        let projection = match projection {
2145            Some(projection) => projection,
2146            None => crate::project_point_to_surface(surface, point).or_refuse(KernelStage::Refine, "project_point_to_surface")?,
2147        };
2148        previous_foot = Some([projection.u, projection.v]);
2149        feet.push([projection.u, projection.v]);
2150        let mut u = projection.u;
2151        if let Some(previous) = previous_u {
2152            while u - previous > period * 0.5 {
2153                u -= period;
2154            }
2155            while previous - u > period * 0.5 {
2156                u += period;
2157            }
2158        }
2159        previous_u = Some(u);
2160        proj_u.push(u);
2161        proj_v.push(projection.v);
2162        point3.push(point);
2163        params.push(fraction);
2164    }
2165    project_timer.stop();
2166    // A biperiodic carrier needs the SAME unwrap in v, or a track that walks up
2167    // to the v seam comes back as 0 at the far end and the fit swings across
2168    // the whole carrier.  Anchor the chain on the first INTERIOR sample and
2169    // pull index 0 onto it afterwards: an endpoint sitting exactly ON the seam
2170    // inverts to either boundary at the solver's whim, and letting that
2171    // coin flip anchor the chain would shift the whole piece a period.
2172    if closed_v {
2173        for index in 2..=dense {
2174            while proj_v[index] - proj_v[index - 1] > period_v * 0.5 {
2175                proj_v[index] -= period_v;
2176            }
2177            while proj_v[index - 1] - proj_v[index] > period_v * 0.5 {
2178                proj_v[index] += period_v;
2179            }
2180        }
2181        while proj_v[0] - proj_v[1] > period_v * 0.5 {
2182            proj_v[0] -= period_v;
2183        }
2184        while proj_v[1] - proj_v[0] > period_v * 0.5 {
2185            proj_v[0] += period_v;
2186        }
2187    }
2188    // Interior mean decides which meridian a crossing endpoint sits on, and
2189    // whether the whole track needs a full-period shift into the domain.
2190    let mean_u: f64 = proj_u[1..dense].iter().sum::<f64>() / (dense - 1) as f64;
2191    let seam_u = if mean_u - u0 > u1 - mean_u { u1 } else { u0 };
2192    // Pin a crossing endpoint onto the seam it actually crosses.  Each
2193    // candidate keeps the OTHER coordinate as projected, so the losing axis
2194    // moves the endpoint bodily off its own 3D vertex — the residual names the
2195    // winner with no surface-type special-casing.  u keeps the whole-track
2196    // mean; a v run can legitimately cover a full period, so its boundary is
2197    // read off the ADJACENT INTERIOR sample instead.  Both pins land in the
2198    // track's own (pre-shift) period, so the shift below carries them along
2199    // with the rest of the track.
2200    let pinned = |index: usize,
2201                  neighbour: usize,
2202                  proj_u: &[f64],
2203                  proj_v: &[f64]|
2204     -> Result<(f64, f64), KernelRefusal> {
2205        let point = point3[index];
2206        let u_error = surface
2207            .evaluate(seam_u, proj_v[index]).or_refuse(KernelStage::Refine, "evaluate")?
2208            .sub(point)
2209            .length();
2210        if closed_v {
2211            let seam_v = nearest_seam_image(proj_v[neighbour], v0, v1);
2212            let v_error = surface
2213                .evaluate(proj_u[index], seam_v).or_refuse(KernelStage::Refine, "evaluate")?
2214                .sub(point)
2215                .length();
2216            if v_error < u_error {
2217                return Ok((proj_u[index], seam_v));
2218            }
2219        }
2220        Ok((seam_u, proj_v[index]))
2221    };
2222    if start_is_crossing {
2223        let (u, v) = pinned(0, 1, &proj_u, &proj_v)?;
2224        proj_u[0] = u;
2225        proj_v[0] = v;
2226    }
2227    if end_is_crossing {
2228        let (u, v) = pinned(dense, dense - 1, &proj_u, &proj_v)?;
2229        proj_u[dense] = u;
2230        proj_v[dense] = v;
2231    }
2232    let mut shift = 0.0;
2233    while mean_u + shift > u1 + 1e-9 {
2234        shift -= period;
2235    }
2236    while mean_u + shift < u0 - 1e-9 {
2237        shift += period;
2238    }
2239    for u in proj_u.iter_mut() {
2240        *u += shift;
2241    }
2242    if closed_v {
2243        let mean_v: f64 = proj_v[1..dense].iter().sum::<f64>() / (dense - 1) as f64;
2244        let mut shift_v = 0.0;
2245        while mean_v + shift_v > v1 + 1e-9 {
2246            shift_v -= period_v;
2247        }
2248        while mean_v + shift_v < v0 - 1e-9 {
2249            shift_v += period_v;
2250        }
2251        for v in proj_v.iter_mut() {
2252            *v += shift_v;
2253        }
2254    }
2255    // Coarsest interpolation whose fitted pcurve reproduces the projected TRACK
2256    // to the kernel's pcurve refinement floor — the floor every other pcurve fit
2257    // is held to, and the one the shell vector-area bar is built from.
2258    //
2259    // It used to stop at 0.002 absolute ("half the validator's tolerance"),
2260    // which is `pcurve_consistency / 2` on any model under 20 units and so four
2261    // decades looser than any other trim: a torus face carrying a collar's rail
2262    // was accepted 1.19e-3 off it, and that one coedge was the whole of a
2263    // 1.94e-4 shell vector-area residual, 14.5x its face's bar.
2264    //
2265    // The deviation is measured against the projected track, not the piece's
2266    // 3D points: a pcurve lies on its surface by construction, so a piece that
2267    // is itself off the carrier would hold the ladder at its top however fine
2268    // the fit, and that is a RAIL defect for the march to answer, not this fit.
2269    let tolerance = crate::pcurve::PCURVE_REFINEMENT_TOLERANCE;
2270    // The ladder is the one verified fitter (`track_fit`), which checks every
2271    // rung at the MIDPOINTS between dense samples, projected on their own —
2272    // checked at the dense samples, the top rung, which interpolates every one
2273    // of them, reads zero by construction and passes whatever lies between (on
2274    // the 574-control-point collar rail it read 9.6e-15) — and breaks the fit,
2275    // C1, where the track crosses a knot line of the carrier.
2276    let at = |fraction: f64| piece.evaluate(d0 + (d1 - d0) * fraction);
2277    let mut track = crate::blend::track_fit::Track {
2278        fractions: params,
2279        uv: proj_u.iter().zip(&proj_v).map(|(u, v)| [*u, *v]).collect(),
2280        feet,
2281        breaks: Vec::new(),
2282    };
2283    let breaks = crate::blend::track_fit::curve_breaks(piece, d0, d1, true);
2284    let polish_timer = crate::blend::BlendTimer::start("timing.track.polish_ms");
2285    crate::blend::track_fit::polish_track(surface, &at, &mut track)?;
2286    polish_timer.stop();
2287    let crossings_timer = crate::blend::BlendTimer::start("timing.track.crossings_ms");
2288    crate::blend::track_fit::insert_knot_crossings(surface, &at, &mut track, &breaks)?;
2289    crossings_timer.stop();
2290    crate::blend::record_count("count.track.breaks", track.breaks.len() as u64);
2291    for &fraction in refine_near {
2292        for _ in 0..2 {
2293            let Some(index) = track.fractions.windows(2).position(|pair| pair[0] <= fraction && fraction <= pair[1]) else { break };
2294            if crate::blend::track_fit::refine_track_intervals(surface, &at, &mut track, &[index], ceiling)?.is_none() {
2295                break;
2296            }
2297        }
2298    }
2299    // Where the top rung still misses the floor, the track is refined ONLY in
2300    // the dense intervals whose out-of-sample checks missed (local insertion
2301    // concentrates the rail's stations, and a uniform dense track can leave one
2302    // sample per knot span where the rail turns hardest), within the dense
2303    // ceiling the track was built under; every other sample, the pinned ends
2304    // and the breaks are kept. The floor and the refusal below are unchanged.
2305    let fit_timer = crate::blend::BlendTimer::start("timing.track.fit_local_ms");
2306    let (fit, inserted) = crate::blend::track_fit::fit_track_locally(
2307        surface,
2308        &mut track,
2309        &at,
2310        tolerance,
2311        REFINE_ROUNDS,
2312        ceiling,
2313    )?;
2314    fit_timer.stop();
2315    crate::blend::record_count("count.track.fit_sections", fit.sections as u64);
2316    crate::blend::record_count("count.track.local_inserted", inserted as u64);
2317    if std::env::var("BREP_DEBUG_TRACK_FIT").is_ok() {
2318        eprintln!("TRACK_FIT closed-chain piece: {inserted} sample(s) inserted locally");
2319        eprintln!(
2320            "TRACK_FIT closed-chain piece: miss {:.3e} at {} samples of {} ({} breaks), on floor {}",
2321            fit.miss,
2322            fit.samples,
2323            track.fractions.len(),
2324            track.breaks.len(),
2325            fit.on_floor
2326        );
2327        // Where the miss sits, against what the piece is made of: an endpoint
2328        // pinned onto a seam reads at fraction 0 or 1; an interior miss names
2329        // the piece's knot span there and how that span compares with the
2330        // piece's mean (local insertion clusters spans at the crotch).
2331        if let Some(worst) = fit.worst {
2332            let parameter = d0 + (d1 - d0) * worst.fraction;
2333            let spans: Vec<[f64; 2]> = piece.knots.windows(2)
2334                .filter(|pair| pair[1] > pair[0] && pair[1] > d0 && pair[0] < d1)
2335                .map(|pair| [pair[0], pair[1]]).collect();
2336            let mean = (d1 - d0) / spans.len().max(1) as f64;
2337            let span = spans.iter().position(|pair| parameter >= pair[0] && parameter <= pair[1]);
2338            let width = span.map(|index| spans[index][1] - spans[index][0]);
2339            eprintln!(
2340                "TRACK_FIT closed-chain worst: fraction {:.9} parameter {:.9} (pinned start {start_is_crossing}, pinned end {end_is_crossing}) \
2341                 knot span {:?} of {} width {:?} vs mean {:.3e}; dense bracket [{:.9}, {:.9}] (interval {})",
2342                worst.fraction, parameter, span, spans.len(), width, mean,
2343                worst.bracket[0], worst.bracket[1], worst.interval
2344            );
2345        }
2346    }
2347    // A fit off the floor is refused, not handed on as the best rung. Before the
2348    // fitter broke at knot lines this ladder returned its best rung silently,
2349    // and on the 2026-09-02 two-torus collar that was 1.07e-7 to 8.46e-7 off at
2350    // the top rung on every collar that builds (measured 2026-09-17); broken,
2351    // each of them reaches the floor. What still misses is a rail the march tore
2352    // (2.7 to 3.3 on the oversized radii, a body refused or rejected anyway).
2353    if !fit.on_floor {
2354        if std::env::var("BREP_DEBUG_TRACK_FIT").is_ok() {
2355            // The piece's two ends against the carrier: a rail end that is
2356            // itself off the carrier holds the ladder however fine the fit.
2357            for (label, t) in [("start", d0), ("end", d1)] {
2358                if let Ok(point) = piece.evaluate(t) {
2359                    let standoff = crate::project_point_to_surface(surface, point).map(|p| p.distance).unwrap_or(f64::NAN);
2360                    eprintln!(
2361                        "TRACK_FIT closed-chain piece {label} ({:.9}, {:.9}, {:.9}) stands {standoff:.3e} off the carrier",
2362                        point.x, point.y, point.z
2363                    );
2364                }
2365            }
2366        }
2367        return Err(KernelRefusal::non_convergence(KernelStage::Refine, CHAIN_PCURVE_FLOOR_WHAT, format!(
2368            "{} a closed chain's support piece misses its pcurve by {:.3e} at {} samples, \
2369             against a floor of {tolerance:.1e}",
2370            crate::blend::PCURVE_OFF_FLOOR,
2371            fit.miss,
2372            fit.samples
2373        )));
2374    }
2375    Ok(fit.curve)
2376}
2377
2378
2379
2380/// Immutable interval input; only coefficients are filled lazily. The Newton
2381/// iterations, foot windows and decline ordering remain in the common reader.
2382pub(in crate::blend) struct PreparedWall<'a> {
2383    fine: &'a [Vec3],
2384    coarse: Vec<Vec3>,
2385    fine_stencils: StencilPreparation,
2386    coarse_stencils: StencilPreparation,
2387}
2388
2389struct StencilPreparation {
2390    finite: bool,
2391    coefficients: std::cell::RefCell<Vec<Option<[Vec3; 4]>>>,
2392}
2393impl StencilPreparation {
2394    fn new(nodes: &[Vec3]) -> Self {
2395        Self { finite: wall_nodes_finite(nodes), coefficients: std::cell::RefCell::new(vec![None; nodes.len().saturating_sub(3)]) }
2396    }
2397    fn coefficients(&self, nodes: &[Vec3], base: usize) -> [Vec3; 4] {
2398        let mut entries = self.coefficients.borrow_mut();
2399        *entries[base].get_or_insert_with(|| wall_coefficients(nodes, base))
2400    }
2401}
2402
2403fn wall_nodes_finite(nodes: &[Vec3]) -> bool {
2404    nodes.iter().all(|c| c.x.is_finite() && c.y.is_finite() && c.z.is_finite())
2405}
2406
2407fn wall_coefficients(nodes: &[Vec3], base: usize) -> [Vec3; 4] {
2408    // Newton's divided differences on 0, 1, 2, 3, s in node units from `base`.
2409    let [p0, p1, p2, p3] = [nodes[base], nodes[base + 1], nodes[base + 2], nodes[base + 3]];
2410    let d1 = p1.sub(p0);
2411    let d2 = p2.sub(p1.scale(2.0)).add(p0).scale(0.5);
2412    let d3 = p3.sub(p2.scale(3.0)).add(p1.scale(3.0)).sub(p0).scale(1.0 / 6.0);
2413    [p0, d1, d2, d3]
2414}
2415
2416fn wall_coarse(fine: &[Vec3]) -> Vec<Vec3> {
2417    fine.iter().step_by(2).copied().collect()
2418}
2419
2420fn wall_preparation_on() -> bool {
2421    std::env::var("BREP_WALL_PREPARE").ok().as_deref() != Some("0")
2422}
2423
2424impl<'a> PreparedWall<'a> {
2425    pub(in crate::blend) fn for_interval(fine: &'a [Vec3]) -> Option<Self> {
2426        (wall_preparation_on() && wall_solve_reuse_on()).then(|| Self::new(fine))
2427    }
2428    fn new(fine: &'a [Vec3]) -> Self {
2429        let coarse = wall_coarse(fine);
2430        let fine_stencils = StencilPreparation::new(fine);
2431        let coarse_stencils = StencilPreparation::new(&coarse);
2432        Self { fine, coarse, fine_stencils, coarse_stencils }
2433    }
2434    pub(in crate::blend) fn read(&self, p: Vec3, j: usize, fraction: f64, radius: f64) -> Option<(f64, f64)> {
2435        let read = |at| declared_wall_probe_prepared(p, self.fine, at, radius, true, Some(self));
2436        match read(j) {
2437            Some(read) => Some(read),
2438            None => match centre_path_node(self.fine, fraction) {
2439                Some(located) if located != j => read(located),
2440                _ => None,
2441            },
2442        }
2443    }
2444}