Skip to main content

brep_kernel/edit/direct_edit/
face_offset_revolution.rs

1use super::*;
2use crate::{ExtendRefusal, KernelRefusal, KernelStage, OrRefuse};
3
4/// Push a general surface of revolution with planar neighbours.
5///
6/// The meridian is densely offset and interpolated, then revolved through the
7/// original full or partial sweep. A dense normal-distance residual check
8/// bounds that approximation. Endpoint rims and start/end meridians are rebuilt,
9/// and planar axial/radial neighbours are re-trimmed around the new boundary.
10/// The offset meridian runs between the source parameters at which it reaches
11/// each rim's cap plane, so a rim is the offset carrier's intersection with its
12/// cap rather than the offset of the old rim (see
13/// [`offset_general_revolution_face`]).
14///
15/// A PARTIAL CYLINDER (an extruded sketch arc, a fillet band between planes)
16/// does not take the general lane: its meridian edges usually meet planes that
17/// are parallel to the axis but do not contain it, where keeping a corner at
18/// its old azimuth is wrong. [`offset_partial_cylinder_face`] re-intersects
19/// those corners exactly first, and only a face outside its scope falls
20/// through to the general lane below.
21///
22/// A boundary whose pcurve is NEITHER a constant-u nor a constant-v line is no
23/// longer refused outright: its 3D image on the offset carrier comes from
24/// [`crate::image_curve`]. That transfer is only valid where the boundary is
25/// definitional, so a general image — or a constant-u meridian — against a
26/// FIXED planar neighbour is measured against that neighbour's plane and
27/// refused, with the number, when it leaves it, because the honest answer
28/// there is a re-intersection this site does not do.
29pub fn offset_revolution_face(
30    solid: &BrepSolid,
31    face_id: u64,
32    distance: f64,
33) -> Result<BrepSolid, KernelRefusal> {
34    if !distance.is_finite() {
35        return Err(KernelRefusal::input(
36            KernelStage::Collect,
37            "distance",
38            "offset_revolution_face: distance must be finite",
39        ));
40    }
41    let scale = solid_model_scale(solid);
42    let tolerance = (scale * 1e-7).max(1e-9);
43    let residual_tolerance = (scale * 1e-4).max(1e-6);
44    // Bar for the general 3D-image ladder: a CURVE-fit accuracy target, which
45    // is a different question from `residual_tolerance`'s bound on the offset
46    // SURFACE, so it takes the kernel's single named fit-accuracy field.
47    let fit_tolerance = crate::KernelTolerances::for_scale(scale, 1e-7).intersection_fit;
48    let (shell_index, face_index) = find_face(solid, face_id).ok_or_else(|| {
49        KernelRefusal::input(
50            KernelStage::Collect,
51            "face_id",
52            format!("offset_revolution_face: no face {face_id}"),
53        )
54    })?;
55    let face = &solid.shells[shell_index].faces[face_index];
56    if !matches!(face.surface.analytic(), Some(AnalyticSurface::Revolution { .. })) {
57        return Err(KernelRefusal::input(
58            KernelStage::Collect,
59            "not_revolution",
60            "offset_revolution_face: the pushed face is not a general revolution",
61        ));
62    }
63    if let Some(result) = offset_partial_cylinder_face(
64        solid,
65        face_id,
66        distance,
67        tolerance,
68        residual_tolerance,
69        scale,
70    )? {
71        return Ok(result);
72    }
73    let trace = BoundaryTrace::open(solid, face_id, distance);
74    let outcome = offset_general_revolution_face(
75        solid,
76        face_id,
77        distance,
78        tolerance,
79        residual_tolerance,
80        fit_tolerance,
81        scale,
82        &trace,
83    );
84    trace.outcome(&outcome);
85    outcome
86}
87
88/// The general lane of [`offset_revolution_face`]: every face the partial
89/// cylinder lane does not take.
90///
91/// Its constant-v RIMS are solved, not transferred. A rim is the circle a
92/// point of the generatrix sweeps, and against a FIXED cap the new rim is the
93/// offset carrier's intersection with that cap's plane. The offset of the
94/// source rim's own generatrix point is not in that plane wherever the
95/// generatrix meets the cap at an angle — the offset normal there has an axial
96/// component `d·n_axis` — and rebuilding the rim at that point left the lib
97/// barrel's rims 8.602e-2 and 6.155e-2 off their caps at d = 0.15 and a solid
98/// that validated 1.757 above the exact offset's volume. So each rim's source
99/// parameter is solved where the offset meridian reaches its cap's plane
100/// ([`solve_rim_span`]), the offset generatrix is sampled between those two
101/// parameters — through the source carrier's natural extension where a solved
102/// parameter leaves its domain — and every rebuilt rim is measured against its
103/// plane and refused, with the number, past the pcurve floor.
104#[allow(clippy::too_many_arguments)]
105fn offset_general_revolution_face(
106    solid: &BrepSolid,
107    face_id: u64,
108    distance: f64,
109    tolerance: f64,
110    residual_tolerance: f64,
111    fit_tolerance: f64,
112    scale: f64,
113    trace: &BoundaryTrace,
114) -> Result<BrepSolid, KernelRefusal> {
115    let (shell_index, face_index) = find_face(solid, face_id).ok_or_else(|| {
116        KernelRefusal::internal(
117            KernelStage::Collect,
118            "face_lookup",
119            format!("offset_revolution_face: no face {face_id}"),
120        )
121    })?;
122    let face = &solid.shells[shell_index].faces[face_index];
123    let Some(AnalyticSurface::Revolution { sweep, .. }) = face.surface.analytic() else {
124        return Err(KernelRefusal::internal(
125            KernelStage::Collect,
126            "carrier_kind",
127            "offset_revolution_face: the pushed face is not a general revolution",
128        ));
129    };
130    let source_structure = crate::revolution_structure(&face.surface).ok_or_else(|| {
131        KernelRefusal::internal(
132            KernelStage::Classify,
133            "revolution_structure",
134            "offset_revolution_face: source carrier lost its revolution structure",
135        )
136    })?;
137    let [u0, u1] = face.surface.domain_u().or_refuse(KernelStage::Classify, "domain_u")?;
138    let [v0, v1] = face.surface.domain_v().or_refuse(KernelStage::Classify, "domain_v")?;
139    if (v0.abs() + (v1 - 1.0).abs()) > 1e-9 {
140        return Err(KernelRefusal::unsupported(
141            KernelStage::Classify,
142            "generatrix_domain",
143            "offset_revolution_face: a non-normalized generatrix domain is deferred",
144        ));
145    }
146    let mut faces_of_edge: HashMap<u64, Vec<u64>> = HashMap::default();
147    for shell in &solid.shells {
148        for candidate in &shell.faces {
149            for loop_record in &candidate.loops {
150                for coedge in &loop_record.coedges {
151                    faces_of_edge
152                        .entry(coedge.edge_id)
153                        .or_default()
154                        .push(candidate.id);
155                }
156            }
157        }
158    }
159    let edge_by_id: HashMap<u64, &EdgeRecord> =
160        solid.edges.iter().map(|edge| (edge.id, edge)).collect();
161    // Offset the meridian densely, then revolve it. A general curve offset is
162    // not exactly rational; using 33 cubic interpolation stations materially
163    // lowers the residual compared with fitting the whole surface back onto its
164    // original sparse control net, while retaining the [0,1] pcurve contract.
165    let u_probe = u0 + (u1 - u0) * 0.137;
166    let span = solve_rim_span(
167        solid,
168        face,
169        face_id,
170        &faces_of_edge,
171        &source_structure.frame,
172        u_probe,
173        distance,
174        residual_tolerance,
175    )?;
176    let source_surface = span.surface.as_ref().unwrap_or(&face.surface);
177    // The dense meridian sample is a pointwise offset — the shared evaluator's
178    // `Face` lane, whose distance is signed ALONG the oriented normal exactly
179    // like direct-edit's own convention (audit §4.1: this side of the kernel
180    // grows the solid on a positive distance).
181    let offsets = OffsetEvaluator::new(
182        "push_revolution",
183        source_surface,
184        OffsetNormal::Face {
185            same_sense: face.same_sense,
186        },
187    );
188    // The carrier's generatrix parameter `w` runs over [0, 1] while the source
189    // parameter it offsets runs over the solved span.
190    let [s0, s1] = span.new;
191    let source_parameter = |w: f64| s0 + (s1 - s0) * w;
192    let mut meridian_points = Vec::with_capacity(33);
193    let mut meridian_parameters = Vec::with_capacity(33);
194    for index in 0..=32 {
195        let parameter = index as f64 / 32.0;
196        let v = source_parameter(v0 + (v1 - v0) * parameter);
197        let target_point = offsets.at(u_probe, v, distance).or_refuse(KernelStage::Refine, "offset_at")?.point;
198        let relative = target_point.sub(source_structure.frame.origin);
199        let axial = relative.dot(source_structure.frame.axis);
200        let radial = relative
201            .sub(source_structure.frame.axis.scale(axial))
202            .length();
203        meridian_points.push(
204            source_structure
205                .frame
206                .origin
207                .add(source_structure.frame.x_axis.scale(radial))
208                .add(source_structure.frame.axis.scale(axial)),
209        );
210        meridian_parameters.push(parameter);
211    }
212    let offset_generatrix = crate::interpolate_curve(&meridian_points, 3, &meridian_parameters).or_refuse(KernelStage::Refine, "interpolate_curve")?;
213    let offset = crate::make_revolution(
214        source_structure.frame.origin,
215        source_structure.frame.axis,
216        &offset_generatrix,
217        *sweep,
218    ).or_refuse(KernelStage::Refine, "make_revolution")?;
219    let structure = crate::revolution_structure(&offset).ok_or_else(|| {
220        KernelRefusal::internal(
221            KernelStage::Refine,
222            "offset_structure",
223            "offset_revolution_face: offset carrier lost its revolution structure — refusing",
224        )
225    })?;
226    if !matches!(offset.analytic(), Some(AnalyticSurface::Revolution { .. })) {
227        return Err(KernelRefusal::internal(
228            KernelStage::Refine,
229            "offset_kind",
230            "offset_revolution_face: offset carrier changed analytic kind — refusing",
231        ));
232    }
233
234    // A general interpolated offset is approximate. Admit it only when its
235    // same-parameter displacement stays normal and at the requested distance.
236    let target = distance.abs();
237    let (mut worst_distance, mut worst_normal) = (0.0f64, 0.0f64);
238    // Stagger away from repeated knots: derivatives at an exact full-period
239    // seam or a multiple internal knot are side-dependent and may be zero in
240    // the generic evaluator even though the carrier is regular there.
241    for iu in 0..13 {
242        for iv in 0..13 {
243            let u = u0 + (u1 - u0) * (iu as f64 + 0.37) / 13.0;
244            let w = v0 + (v1 - v0) * (iv as f64 + 0.41) / 13.0;
245            let v = source_parameter(w);
246            let source_point = source_surface.evaluate(u, v).or_refuse(KernelStage::Refine, "evaluate")?;
247            let offset_point = offset.evaluate(u, w).or_refuse(KernelStage::Refine, "evaluate")?;
248            let delta = offset_point.sub(source_point);
249            let normal = offsets
250                .normal(u, v)
251                .map_err(|error| {
252                    KernelRefusal::internal(
253                        KernelStage::Refine,
254                        "sample_normal",
255                        format!("offset_revolution_face: sample normal: {error}"),
256                    )
257                })?;
258            let distance_error = (delta.length() - target).abs();
259            let normal_error = if delta.length() > tolerance {
260                1.0 - delta
261                    .normalized()
262                    .map_err(|error| {
263                        KernelRefusal::internal(
264                            KernelStage::Refine,
265                            "sample_delta",
266                            format!("offset_revolution_face: sample delta: {error}"),
267                        )
268                    })?
269                    .dot(normal)
270                    .abs()
271            } else {
272                0.0
273            };
274            if distance_error > residual_tolerance || normal_error > 5e-4 {
275                return Err(KernelRefusal::non_convergence(
276                    KernelStage::Refine,
277                    "offset_fit",
278                    format!(
279                    "offset_revolution_face: offset fit exceeds tolerance \
280                     (distance error {distance_error:.3e}, normal error {normal_error:.3e})"
281                ),
282                ));
283            }
284            worst_distance = worst_distance.max(distance_error);
285            worst_normal = worst_normal.max(normal_error);
286        }
287    }
288    trace.write(format!(
289        "carrier\trims_old={:?}\trims_new={:?}\textended={}\t\
290         fit_distance_error={worst_distance:.3e}\tfit_normal_error={worst_normal:.3e}",
291        span.old,
292        span.new,
293        span.surface.is_some()
294    ));
295
296    // The offset generatrix endpoints revolve into the only two possible rims:
297    // the low rim at w = 0, the high rim at w = 1.
298    let [g0, g1] = structure.generatrix.domain().or_refuse(KernelStage::Refine, "domain")?;
299    let mut rim_candidates = Vec::new();
300    for parameter in [g0, g1] {
301        let point = structure.generatrix.evaluate(parameter).or_refuse(KernelStage::Refine, "evaluate")?;
302        let relative = point.sub(structure.frame.origin);
303        let axial = relative.dot(structure.frame.axis);
304        let radial = relative.sub(structure.frame.axis.scale(axial));
305        let radius = radial.length();
306        if radius <= tolerance {
307            return Err(KernelRefusal::input(
308                KernelStage::Refine,
309                "rim_collapse",
310                format!(
311                "offset_revolution_face: an offset rim collapses onto the axis \
312                 (point {point:?}, frame {:?}) — refusing",
313                structure.frame
314            ),
315            ));
316        }
317        rim_candidates.push(crate::make_arc(
318            structure
319                .frame
320                .origin
321                .add(structure.frame.axis.scale(axial)),
322            structure.frame.x_axis,
323            structure.frame.y_axis,
324            radius,
325            0.0,
326            structure.sweep,
327        ).or_refuse(KernelStage::Refine, "make_arc")?);
328    }
329
330    // Each boundary's pcurve on the new carrier. A rim that is not at a
331    // generatrix end — a face trimmed across its carrier — moves to w = 0 or 1,
332    // so the ISO boundaries are carried by the affine map that takes the two
333    // rims' old parameters there. A GENERAL boundary is a transfer of the old
334    // one, so it keeps its source parameter: w = (s − s₀)/(s₁ − s₀).
335    let uv_tolerance = 1e-7;
336    let pcurves = face
337        .loops
338        .iter()
339        .map(|loop_record| {
340            loop_record
341                .coedges
342                .iter()
343                .map(|coedge| match iso_kind(&coedge.pcurve, uv_tolerance)? {
344                    Some(_) => span.carry_pcurve(&coedge.pcurve),
345                    None => span.transfer_pcurve(&coedge.pcurve),
346                })
347                .collect::<Result<Vec<_>, KernelRefusal>>()
348        })
349        .collect::<Result<Vec<_>, KernelRefusal>>()?;
350
351    let floor = crate::pcurve::PCURVE_REFINEMENT_TOLERANCE;
352    let mut new_curves: HashMap<u64, NurbsCurve> = HashMap::default();
353    let mut new_vertices: HashMap<u64, Vec3> = HashMap::default();
354    let mut planar_neighbours: HashSet<u64> = HashSet::default();
355    for (loop_record, loop_pcurves) in face.loops.iter().zip(&pcurves) {
356        for (coedge, pcurve) in loop_record.coedges.iter().zip(loop_pcurves) {
357            if new_curves.contains_key(&coedge.edge_id) {
358                continue;
359            }
360            let edge = *edge_by_id.get(&coedge.edge_id).ok_or_else(|| {
361                KernelRefusal::internal(
362                    KernelStage::Classify,
363                    "missing_edge",
364                    format!("offset_revolution_face: missing edge {}", coedge.edge_id),
365                )
366            })?;
367            let incident = faces_of_edge.get(&edge.id).cloned().unwrap_or_default();
368            let neighbour = incident.iter().find(|candidate| **candidate != face_id);
369            // Kept for the gates below: an edge shared with a FIXED planar
370            // cap has to stay in that cap's plane, and only this loop knows
371            // which plane that is.
372            let mut neighbour_plane = None;
373            if let Some(neighbour) = neighbour {
374                let (ns, nf) = find_face(solid, *neighbour).ok_or_else(|| {
375                    KernelRefusal::internal(
376                        KernelStage::Classify,
377                        "missing_neighbour",
378                        format!("offset_revolution_face: missing neighbour {neighbour}"),
379                    )
380                })?;
381                match plane_of_surface(
382                    &solid.shells[ns].faces[nf].surface,
383                    residual_tolerance,
384                    "offset_revolution_face",
385                ) {
386                    Err(_) => {
387                        return Err(KernelRefusal::unsupported(
388                            KernelStage::Classify,
389                            "curved_neighbour",
390                            "offset_revolution_face: only planar neighbours are supported \
391                                for a general revolution — refusing",
392                        ));
393                    }
394                    Ok(plane) => neighbour_plane = Some(plane),
395                }
396                planar_neighbours.insert(*neighbour);
397            }
398
399            // Constant-v pcurve → a revolved endpoint rim. Constant-u → a
400            // start/end meridian (or the intrinsic full-sweep seam).
401            let mut general_image = false;
402            let mut meridian = false;
403            let mut curve = match iso_kind(&coedge.pcurve, uv_tolerance)? {
404                Some(IsoKind::ConstantV) => {
405                    let [q0, _] = pcurve.domain().or_refuse(KernelStage::Refine, "domain")?;
406                    let w = pcurve.evaluate(q0).or_refuse(KernelStage::Refine, "evaluate")?.y;
407                    if w.abs() <= uv_tolerance {
408                        rim_candidates[0].clone()
409                    } else if (w - 1.0).abs() <= uv_tolerance {
410                        rim_candidates[1].clone()
411                    } else {
412                        return Err(KernelRefusal::unsupported(
413                            KernelStage::Classify,
414                            "interior_rim",
415                            format!(
416                            "offset_revolution_face: constant-v boundary edge {} sits at v = {w} \
417                             between the face's two rims — refusing",
418                            edge.id
419                        ),
420                        ));
421                    }
422                }
423                Some(IsoKind::ConstantU(u)) => {
424                    meridian = true;
425                    offset.iso_curve_u(u).or_refuse(KernelStage::Refine, "iso_curve_u")?
426                }
427                None => {
428                    // A GENERAL boundary pcurve — what this site used to refuse.
429                    // Its image on the offset carrier is built by the shared
430                    // ladder, oriented off `coedge.forward` (the authoritative
431                    // pcurve↔edge relation) rather than the iso lanes' geometric
432                    // nearest-endpoint match, because a general image need not
433                    // share an endpoint with the old curve at all.
434                    general_image = true;
435                    let image = crate::image_curve(
436                        &offset,
437                        pcurve,
438                        fit_tolerance,
439                        "offset_revolution_face",
440                    ).or_refuse(KernelStage::Refine, "image_curve")?;
441                    if coedge.forward {
442                        image.curve
443                    } else {
444                        image.curve.reversed().or_refuse(KernelStage::Refine, "reversed")?
445                    }
446                }
447            };
448            let kind = if general_image {
449                "general"
450            } else if meridian {
451                "meridian"
452            } else {
453                "rim"
454            };
455            trace.edge(kind, edge, neighbour, neighbour_plane.as_ref(), &curve);
456            // Every rebuilt boundary is measured against the FIXED planar
457            // neighbour it must lie in; the validator's `pcurve_acceptance`
458            // band (2.5% of the model diagonal) is far too loose to catch one
459            // that leaves it.
460            //
461            // A RIM is built in its cap's plane by `solve_rim_span`, so past
462            // the pcurve floor it is a cap this lane cannot meet — one not
463            // normal to the axis, or two coplanar caps of one rim that
464            // disagree — and it is refused with the number.
465            //
466            // A general image, or a constant-u MERIDIAN, is a transfer, which
467            // is only legitimate where the boundary is DEFINITIONAL. Keeping a
468            // meridian's azimuth is exact only when the neighbour plane contains
469            // the axis (a partial revolve's radial end face); a plane beside the
470            // axis — the straight sketch segment next to an extruded arc —
471            // needs the corner re-intersected, and the kept azimuth left it
472            // 0.33 off that plane on a reported extrude, where re-trimming the
473            // neighbour around it tilted the wall and still validated.
474            if let Some(plane) = &neighbour_plane {
475                let [g0, g1] = curve.domain().or_refuse(KernelStage::Refine, "domain")?;
476                let mut worst = 0.0f64;
477                for index in 0..=32 {
478                    let t = g0 + (g1 - g0) * index as f64 / 32.0;
479                    let point = curve.evaluate(t).or_refuse(KernelStage::Refine, "evaluate")?;
480                    worst = worst.max(point.sub(plane.origin).dot(plane.normal).abs());
481                }
482                if !(general_image || meridian) && !(worst <= floor) {
483                    return Err(KernelRefusal::unsupported(
484                        KernelStage::Refine,
485                        "rim_off_cap",
486                        format!(
487                        "offset_revolution_face: the rebuilt rim of boundary edge {} leaves the \
488                         plane of neighbour face {} by {worst:.3e} (pcurve floor {floor:.0e}) — \
489                         refusing",
490                        edge.id,
491                        neighbour.copied().unwrap_or_default()
492                    ),
493                    ));
494                }
495                if (general_image || meridian) && worst > residual_tolerance {
496                    let lane = if meridian {
497                        "meridian at its old azimuth"
498                    } else {
499                        "general image"
500                    };
501                    return Err(KernelRefusal::unsupported(
502                        KernelStage::Refine,
503                        "transfer_off_neighbour",
504                        format!(
505                        "offset_revolution_face: the {lane} of boundary edge {} \
506                         leaves its fixed planar neighbour by {worst:.3e} (limit \
507                         {residual_tolerance:.3e}) — re-intersecting the offset carrier \
508                         with a fixed neighbour is a separate capability; refusing",
509                        edge.id
510                    ),
511                    ));
512                }
513            }
514            if !general_image {
515                let [c0, c1] = curve.domain().or_refuse(KernelStage::Refine, "domain")?;
516                let old_start = edge.curve.evaluate(edge.t0).or_refuse(KernelStage::Refine, "evaluate")?;
517                if curve.evaluate(c0).or_refuse(KernelStage::Refine, "evaluate")?.sub(old_start).length()
518                    > curve.evaluate(c1).or_refuse(KernelStage::Refine, "evaluate")?.sub(old_start).length()
519                {
520                    curve = curve.reversed().or_refuse(KernelStage::Refine, "reversed")?;
521                }
522            }
523            // A corner is rebuilt by each boundary edge that ends on it, and
524            // those must agree: a rim and a meridian meeting at one vertex, or
525            // two rim arcs of one split circle. One that does not is a rim this
526            // lane built across a vertex it does not end at.
527            let [d0, d1] = curve.domain().or_refuse(KernelStage::Refine, "domain")?;
528            let ends = [
529                (edge.start_vertex_id, curve.evaluate(d0).or_refuse(KernelStage::Refine, "evaluate")?),
530                (edge.end_vertex_id, curve.evaluate(d1).or_refuse(KernelStage::Refine, "evaluate")?),
531            ];
532            for (vertex_id, point) in ends {
533                if let Some(previous) = new_vertices.insert(vertex_id, point) {
534                    let gap = previous.sub(point).length();
535                    if !(gap <= floor) {
536                        return Err(KernelRefusal::internal(
537                            KernelStage::Fragment,
538                            "vertex_agreement",
539                            format!(
540                            "offset_revolution_face: boundary edge {} rebuilds vertex {vertex_id} \
541                             {gap:.3e} from where another boundary edge put it (pcurve floor \
542                             {floor:.0e}) — refusing",
543                            edge.id
544                        ),
545                        ));
546                    }
547                }
548            }
549            new_curves.insert(edge.id, curve);
550        }
551    }
552
553    let mut result = solid.clone();
554    let pushed = &mut result.shells[shell_index].faces[face_index];
555    if span.old != [0.0, 1.0] || span.new != [0.0, 1.0] {
556        for (loop_record, loop_pcurves) in pushed.loops.iter_mut().zip(pcurves) {
557            for (coedge, pcurve) in loop_record.coedges.iter_mut().zip(loop_pcurves) {
558                coedge.pcurve = pcurve;
559            }
560        }
561    }
562    pushed.surface = offset;
563    retrim_planar_neighbours(
564        solid,
565        result,
566        &edge_by_id,
567        &new_curves,
568        &new_vertices,
569        planar_neighbours,
570        residual_tolerance,
571        scale,
572    )
573}
574
575/// Where the offset carrier's generatrix starts and ends.
576struct RimSpan {
577    /// The face's low and high rim parameters on the SOURCE generatrix — the
578    /// domain end where the face has no rim on that side.
579    old: [f64; 2],
580    /// The source parameters whose offset meridian points lie in those rims'
581    /// cap planes: the new carrier's w = 0 and w = 1.
582    new: [f64; 2],
583    /// The source carrier extended past its domain by its own terminal span
584    /// ([`NurbsSurface::extend_natural`]), when a solved parameter leaves it.
585    surface: Option<NurbsSurface>,
586}
587
588impl RimSpan {
589    /// An iso boundary's pcurve on the new carrier: `v ↦ (v − old₀)/(old₁ −
590    /// old₀)`, which puts the two rims at w = 0 and 1.
591    fn carry_pcurve(&self, pcurve: &NurbsCurve) -> Result<NurbsCurve, KernelRefusal> {
592        Self::map_v(pcurve, self.old)
593    }
594
595    /// A general boundary's pcurve on the new carrier: `v ↦ (v − new₀)/(new₁ −
596    /// new₀)`, the same source parameter and so the same offset point.
597    fn transfer_pcurve(&self, pcurve: &NurbsCurve) -> Result<NurbsCurve, KernelRefusal> {
598        Self::map_v(pcurve, self.new)
599    }
600
601    /// `pcurve` with `v ↦ (v − a₀)/(a₁ − a₀)` applied to its control net (an
602    /// affine map, so exact), or the pcurve itself when that is the identity.
603    fn map_v(pcurve: &NurbsCurve, [a0, a1]: [f64; 2]) -> Result<NurbsCurve, KernelRefusal> {
604        if [a0, a1] == [0.0, 1.0] {
605            return Ok(pcurve.clone());
606        }
607        let mut carried = pcurve.clone();
608        for control in &mut carried.control_points {
609            let point = control.point().or_refuse(KernelStage::Classify, "point")?;
610            *control = crate::Vec4::from_point(
611                Vec3::new(point.x, (point.y - a0) / (a1 - a0), point.z),
612                control.w,
613            );
614        }
615        Ok(carried)
616    }
617}
618
619/// Solve the source generatrix parameters between which the offset carrier
620/// must run so that each constant-v rim of `face` lies in the plane of the
621/// fixed cap across it.
622///
623/// A rim at source parameter `a` against a cap whose plane the axis pierces at
624/// station `h` moves to the root of `f(s) = axial(S(u, s) + d·N(u, s)) − h`
625/// nearest `a`. The root is bracketed by stepping from `a` along the source
626/// meridian's own axial rate (doubling the step up to four times), extending
627/// the source carrier by its natural continuation where the bracket leaves the
628/// domain, and bisected to exhaustion. A rim whose generatrix leaves its cap
629/// tangentially (no axial rate to step along, and a crossing on either side of
630/// the rim), and one whose root cannot be bracketed, are refused by name with
631/// the drift.
632#[allow(clippy::too_many_arguments)]
633fn solve_rim_span(
634    solid: &BrepSolid,
635    face: &FaceRecord,
636    face_id: u64,
637    faces_of_edge: &HashMap<u64, Vec<u64>>,
638    frame: &crate::RevolutionFrame,
639    u: f64,
640    distance: f64,
641    residual_tolerance: f64,
642) -> Result<RimSpan, KernelRefusal> {
643    let op = "offset_revolution_face";
644    let uv_tolerance = 1e-7;
645    // (source parameter, cap station, cap face) of every constant-v boundary
646    // with a planar neighbour, and the face's parameter extent.
647    let mut rims: Vec<(f64, f64, u64)> = Vec::new();
648    let (mut v_min, mut v_max) = (f64::INFINITY, f64::NEG_INFINITY);
649    for coedge in face.loops.iter().flat_map(|loop_record| &loop_record.coedges) {
650        let [q0, q1] = coedge.pcurve.domain().or_refuse(KernelStage::Classify, "domain")?;
651        let uv0 = coedge.pcurve.evaluate(q0).or_refuse(KernelStage::Classify, "evaluate")?;
652        let uv1 = coedge.pcurve.evaluate(q1).or_refuse(KernelStage::Classify, "evaluate")?;
653        let uvm = coedge.pcurve.evaluate(0.5 * (q0 + q1)).or_refuse(KernelStage::Classify, "evaluate")?;
654        for uv in [uv0, uv1, uvm] {
655            v_min = v_min.min(uv.y);
656            v_max = v_max.max(uv.y);
657        }
658        if !matches!(iso_kind(&coedge.pcurve, uv_tolerance)?, Some(IsoKind::ConstantV)) {
659            continue;
660        }
661        let incident = faces_of_edge.get(&coedge.edge_id).cloned().unwrap_or_default();
662        let Some(neighbour) = incident.iter().find(|candidate| **candidate != face_id) else {
663            continue;
664        };
665        let (ns, nf) = find_face(solid, *neighbour).ok_or_else(|| {
666            KernelRefusal::internal(
667                KernelStage::Classify,
668                "rim_neighbour",
669                format!("{op}: missing neighbour {neighbour}"),
670            )
671        })?;
672        // A curved cap is refused by the boundary loop, ahead of any rim.
673        let Ok(plane) =
674            plane_of_surface(&solid.shells[ns].faces[nf].surface, residual_tolerance, op)
675        else {
676            continue;
677        };
678        // Where the axis pierces the cap's plane.
679        let station =
680            plane.origin.sub(frame.origin).dot(plane.normal) / frame.axis.dot(plane.normal);
681        rims.push((uv0.y, station, *neighbour));
682    }
683    if rims.is_empty() {
684        return Ok(RimSpan {
685            old: [0.0, 1.0],
686            new: [0.0, 1.0],
687            surface: None,
688        });
689    }
690    // One rim per end: every constant-v boundary sits at the face's lowest or
691    // highest parameter.
692    let mut ends: [Option<(f64, f64, u64)>; 2] = [None, None];
693    for rim in &rims {
694        let slot = if (rim.0 - v_min).abs() <= (rim.0 - v_max).abs() { 0 } else { 1 };
695        match ends[slot] {
696            None => ends[slot] = Some(*rim),
697            Some(first) if (first.0 - rim.0).abs() > uv_tolerance => {
698                return Err(KernelRefusal::unsupported(
699                    KernelStage::Classify,
700                    "two_rims_one_side",
701                    format!(
702                    "{op}: face {face_id} has constant-v boundaries at v = {} and v = {} on one \
703                     side of its trim — refusing",
704                    first.0, rim.0
705                ),
706                ));
707            }
708            Some(_) => {}
709        }
710    }
711
712    let same_sense = face.same_sense;
713    let mut surface: Option<NurbsSurface> = None;
714    let mut old = [0.0, 1.0];
715    let mut new = [0.0, 1.0];
716    for (slot, end) in ends.iter().enumerate() {
717        let Some((start, station, neighbour)) = *end else {
718            continue;
719        };
720        old[slot] = start;
721        let axial_miss = |carrier: &NurbsSurface, s: f64| -> Result<f64, KernelRefusal> {
722            let convention = OffsetNormal::Face { same_sense };
723            let point = OffsetEvaluator::new("push_revolution", carrier, convention)
724                .at(u, s, distance).or_refuse(KernelStage::Refine, "offset_at")?
725                .point;
726            Ok(point.sub(frame.origin).dot(frame.axis) - station)
727        };
728        let base = surface.clone().unwrap_or_else(|| face.surface.clone());
729        let drift = axial_miss(&base, start)?;
730        if drift == 0.0 {
731            new[slot] = start;
732            continue;
733        }
734        let rate = face.surface.derivatives(u, start, 1).or_refuse(KernelStage::Refine, "derivatives")?[0][1].dot(frame.axis);
735        let step = -drift / rate;
736        if !step.is_finite() {
737            // The generatrix leaves its cap along the cap: the offset meridian
738            // dips through the plane on BOTH sides of the rim, and which
739            // crossing is the rim is not the solve's to pick.
740            return Err(KernelRefusal::ill_posed(
741                KernelStage::Refine,
742                "tangent_rim",
743                format!(
744                "{op}: the generatrix of face {face_id} meets the plane of neighbour face \
745                 {neighbour} tangentially at its rim, which the push leaves {drift:.3e} off \
746                 that plane — refusing"
747            ),
748            ));
749        }
750        let mut bracket = None;
751        let mut reach = 2.0 * step;
752        for _ in 0..4 {
753            if !reach.is_finite() {
754                break;
755            }
756            let far = start + reach;
757            let carrier = extend_to(&base, far)
758                .map_err(|refusal| {
759                    // The extension's own text rides along verbatim; its class is read off
760                    // the variant — a lower-level failure is a defect, every other variant
761                    // is a limit of the natural continuation this lane cannot cross.
762                    let message = {
763                        let refusal = refusal.describe();
764                        format!(
765                        "{op}: the rim of face {face_id} against neighbour face {neighbour} lies \
766                         {drift:.3e} off that cap, and the carrier cannot be extended to reach \
767                         it: {refusal} — refusing"
768                    )
769                    };
770                    match refusal {
771                        ExtendRefusal::Failed(_) => KernelRefusal::internal(
772                            KernelStage::Refine,
773                            "carrier_extension_failed",
774                            message,
775                        ),
776                        _ => KernelRefusal::unsupported(
777                            KernelStage::Refine,
778                            "carrier_extension",
779                            message,
780                        ),
781                    }
782                })?;
783            let far_miss = axial_miss(carrier.as_ref().unwrap_or(&base), far)?;
784            if far_miss == 0.0 || far_miss.signum() != drift.signum() {
785                bracket = Some((far, carrier));
786                break;
787            }
788            reach *= 2.0;
789        }
790        let Some((far, carrier)) = bracket else {
791            return Err(KernelRefusal::non_convergence(
792                KernelStage::Refine,
793                "rim_bracket",
794                format!(
795                "{op}: the offset meridian of face {face_id} does not reach the plane of \
796                 neighbour face {neighbour} near its rim, which the push leaves {drift:.3e} off \
797                 that plane — refusing"
798            ),
799            ));
800        };
801        let carrier_ref = carrier.as_ref().unwrap_or(&base);
802        let (mut near, mut far) = (start, far);
803        for _ in 0..200 {
804            let middle = 0.5 * (near + far);
805            if middle == near || middle == far {
806                break;
807            }
808            if axial_miss(carrier_ref, middle)?.signum() == drift.signum() {
809                near = middle;
810            } else {
811                far = middle;
812            }
813        }
814        new[slot] = 0.5 * (near + far);
815        if carrier.is_some() {
816            surface = carrier;
817        }
818    }
819    if !(new[1] > new[0]) {
820        return Err(KernelRefusal::input(
821            KernelStage::Refine,
822            "rims_crossed",
823            format!(
824            "{op}: pushing face {face_id} by {distance} crosses its two rims over (source \
825             parameters {} and {}) — refusing",
826            new[0], new[1]
827        ),
828        ));
829    }
830    Ok(RimSpan { old, new, surface })
831}
832
833/// Which iso line of its carrier a boundary pcurve is, by the same three-sample
834/// test the lane has always used; `None` for a general pcurve.
835enum IsoKind {
836    ConstantV,
837    ConstantU(f64),
838}
839
840fn iso_kind(pcurve: &NurbsCurve, uv_tolerance: f64) -> Result<Option<IsoKind>, KernelRefusal> {
841    let [q0, q1] = pcurve.domain().or_refuse(KernelStage::Classify, "domain")?;
842    let uv0 = pcurve.evaluate(q0).or_refuse(KernelStage::Classify, "evaluate")?;
843    let uv1 = pcurve.evaluate(q1).or_refuse(KernelStage::Classify, "evaluate")?;
844    let uvm = pcurve.evaluate(0.5 * (q0 + q1)).or_refuse(KernelStage::Classify, "evaluate")?;
845    Ok(
846        if (uv0.y - uv1.y).abs() <= uv_tolerance && (uv0.y - uvm.y).abs() <= uv_tolerance {
847            Some(IsoKind::ConstantV)
848        } else if (uv0.x - uv1.x).abs() <= uv_tolerance && (uv0.x - uvm.x).abs() <= uv_tolerance {
849            Some(IsoKind::ConstantU(uv0.x))
850        } else {
851            None
852        },
853    )
854}
855
856/// `surface` naturally extended in v so that its domain contains `parameter`,
857/// or `None` when it already does.
858fn extend_to(surface: &NurbsSurface, parameter: f64) -> Result<Option<NurbsSurface>, ExtendRefusal> {
859    let [v0, v1] = surface.domain_v().map_err(ExtendRefusal::Failed)?;
860    let (side, delta) = if parameter < v0 {
861        (crate::SurfaceSide::VMin, v0 - parameter)
862    } else if parameter > v1 {
863        (crate::SurfaceSide::VMax, parameter - v1)
864    } else {
865        return Ok(None);
866    };
867    // The solve bisects inside the extension, so its length only has to
868    // clear the knot identity bar `extend_natural` refuses an increment under.
869    surface
870        .extend_natural(side, delta.max(2.0 * crate::curve::KNOT_IDENTITY_TOL))
871        .map(Some)
872}
873
874/// A meridian boundary of a partial cylinder and the fixed plane across it.
875struct Meridian<'a> {
876    edge: &'a EdgeRecord,
877    neighbour: u64,
878    /// Azimuth of the meridian in the carrier's frame.
879    azimuth: f64,
880    plane_origin: Vec3,
881    plane_normal: Vec3,
882}
883
884/// Push a PARTIAL CYLINDER exactly: a straight generatrix parallel to the axis,
885/// swept through less than a full turn, bounded by two meridians and two rims,
886/// with a planar neighbour across each.
887///
888/// The offset carrier is the coaxial cylinder of radius ρ′ = ρ ± d, and every
889/// boundary is re-derived as that carrier's intersection with the fixed
890/// neighbour's carrier:
891/// * a rim against an axis-normal cap stays at its axial station;
892/// * a meridian against a plane PARALLEL to the axis moves to where the new
893///   circle crosses that plane, `ρ′·A·cos(θ − φ) = D` in the carrier frame. A
894///   radial end face (D = 0) keeps its azimuth — the one case the general
895///   lane's kept azimuth is right for.
896///
897/// Of the two crossings, the corner is the one from which the new arc leaves
898/// into the side of the plane the old arc left into: the old arc's tangent side
899/// at a transversal corner, the centre's side at a tangent one (a fillet band,
900/// a sketch arc tangent to its segment), where "nearest crossing" is a tie.
901///
902/// `Ok(None)` means the face is not this lane's shape and goes to the general
903/// lane; an `Err` is a push this lane recognised and refuses by name.
904fn offset_partial_cylinder_face(
905    solid: &BrepSolid,
906    face_id: u64,
907    distance: f64,
908    tolerance: f64,
909    residual_tolerance: f64,
910    scale: f64,
911) -> Result<Option<BrepSolid>, KernelRefusal> {
912    use std::f64::consts::TAU;
913    let op = "offset_revolution_face";
914    let (shell_index, face_index) = find_face(solid, face_id).ok_or_else(|| {
915        KernelRefusal::internal(
916            KernelStage::Collect,
917            "partial_face_lookup",
918            format!("{op}: no face {face_id}"),
919        )
920    })?;
921    let face = &solid.shells[shell_index].faces[face_index];
922    let Some(AnalyticSurface::Revolution {
923        frame, generatrix, ..
924    }) = face.surface.analytic()
925    else {
926        return Ok(None);
927    };
928    if generatrix.degree != 1 || generatrix.control_points.len() != 2 {
929        return Ok(None);
930    }
931    let (origin, axis, x_axis, y_axis) = (frame.origin, frame.axis, frame.x_axis, frame.y_axis);
932    // (axial coordinate, radial vector) of a point about the carrier axis.
933    let cylindrical = |point: Vec3| {
934        let relative = point.sub(origin);
935        let axial = relative.dot(axis);
936        (axial, relative.sub(axis.scale(axial)))
937    };
938    let azimuth_of = |radial: Vec3| radial.dot(y_axis).atan2(radial.dot(x_axis));
939    let (generatrix_z0, generatrix_radial0) = cylindrical(generatrix.control_points[0].point().or_refuse(KernelStage::Classify, "point")?);
940    let (generatrix_z1, generatrix_radial1) = cylindrical(generatrix.control_points[1].point().or_refuse(KernelStage::Classify, "point")?);
941    if generatrix_radial0.sub(generatrix_radial1).length() > residual_tolerance {
942        return Ok(None); // a cone
943    }
944    let radius = generatrix_radial0.length();
945    let span = (generatrix_z1 - generatrix_z0).abs();
946
947    if face.loops.len() != 1 || face.loops[0].coedges.len() != 4 {
948        return Ok(None);
949    }
950    let mut faces_of_edge: HashMap<u64, Vec<u64>> = HashMap::default();
951    for shell in &solid.shells {
952        for candidate in &shell.faces {
953            for loop_record in &candidate.loops {
954                for coedge in &loop_record.coedges {
955                    faces_of_edge
956                        .entry(coedge.edge_id)
957                        .or_default()
958                        .push(candidate.id);
959                }
960            }
961        }
962    }
963    let edge_by_id: HashMap<u64, &EdgeRecord> =
964        solid.edges.iter().map(|edge| (edge.id, edge)).collect();
965
966    // Classify the four boundary edges by their 3D geometry.
967    let mut meridians: Vec<Meridian> = Vec::new();
968    let mut rims: Vec<(&EdgeRecord, f64)> = Vec::new();
969    let mut planar_neighbours: HashSet<u64> = HashSet::default();
970    for coedge in &face.loops[0].coedges {
971        let edge = *edge_by_id
972            .get(&coedge.edge_id)
973            .ok_or_else(|| {
974                KernelRefusal::internal(
975                    KernelStage::Classify,
976                    "partial_missing_edge",
977                    format!("{op}: missing edge {}", coedge.edge_id),
978                )
979            })?;
980        let incident = faces_of_edge.get(&edge.id).cloned().unwrap_or_default();
981        if incident.len() != 2 || incident.iter().filter(|f| **f == face_id).count() != 1 {
982            return Ok(None);
983        }
984        let neighbour = *incident.iter().find(|f| **f != face_id).expect("two incidences");
985        let (ns, nf) = find_face(solid, neighbour).ok_or_else(|| {
986            KernelRefusal::internal(
987                KernelStage::Classify,
988                "partial_missing_neighbour",
989                format!("{op}: missing neighbour {neighbour}"),
990            )
991        })?;
992        let Ok(plane) =
993            plane_of_surface(&solid.shells[ns].faces[nf].surface, residual_tolerance, op)
994        else {
995            return Ok(None);
996        };
997        let (z_start, radial_start) = cylindrical(edge.curve.evaluate(edge.t0).or_refuse(KernelStage::Classify, "evaluate")?);
998        let (z_mid, radial_mid) = cylindrical(edge.curve.evaluate(0.5 * (edge.t0 + edge.t1)).or_refuse(KernelStage::Classify, "evaluate")?);
999        let (z_end, radial_end) = cylindrical(edge.curve.evaluate(edge.t1).or_refuse(KernelStage::Classify, "evaluate")?);
1000        if radial_start.sub(radial_mid).length() <= residual_tolerance
1001            && radial_start.sub(radial_end).length() <= residual_tolerance
1002        {
1003            if plane.normal.dot(axis).abs() * span > residual_tolerance {
1004                return Ok(None); // a plane that is not parallel to the axis
1005            }
1006            meridians.push(Meridian {
1007                edge,
1008                neighbour,
1009                azimuth: azimuth_of(radial_start),
1010                plane_origin: plane.origin,
1011                plane_normal: plane.normal,
1012            });
1013        } else if (z_start - z_mid).abs() <= residual_tolerance
1014            && (z_start - z_end).abs() <= residual_tolerance
1015        {
1016            if plane.normal.cross(axis).length() * radius > residual_tolerance {
1017                return Ok(None); // a cap that is not normal to the axis
1018            }
1019            rims.push((edge, z_start));
1020        } else {
1021            return Ok(None);
1022        }
1023        planar_neighbours.insert(neighbour);
1024    }
1025    if meridians.len() != 2 || rims.len() != 2 || (rims[0].1 - rims[1].1).abs() <= residual_tolerance
1026    {
1027        return Ok(None);
1028    }
1029    // Every corner is a meridian end at one of the two rim stations, each rim
1030    // runs between two corners, and every other edge ending at a corner is a
1031    // straight edge between two of the planes this push re-trims — so the chord
1032    // the re-trim rebuilds stays on both.
1033    let mut stations: HashMap<u64, f64> = HashMap::default();
1034    for meridian in &meridians {
1035        for vertex_id in [meridian.edge.start_vertex_id, meridian.edge.end_vertex_id] {
1036            let vertex = solid
1037                .vertices
1038                .iter()
1039                .find(|vertex| vertex.id == vertex_id)
1040                .ok_or_else(|| {
1041                    KernelRefusal::internal(
1042                        KernelStage::Classify,
1043                        "partial_missing_vertex",
1044                        format!("{op}: missing vertex {vertex_id}"),
1045                    )
1046                })?;
1047            let (axial, _) = cylindrical(vertex.point);
1048            let Some(station) = rims
1049                .iter()
1050                .map(|(_, station)| *station)
1051                .find(|station| (axial - station).abs() <= residual_tolerance)
1052            else {
1053                return Ok(None);
1054            };
1055            stations.insert(vertex_id, station);
1056        }
1057    }
1058    let boundary: HashSet<u64> =
1059        face.loops[0].coedges.iter().map(|coedge| coedge.edge_id).collect();
1060    for edge in &solid.edges {
1061        let at_corner = stations.contains_key(&edge.start_vertex_id)
1062            || stations.contains_key(&edge.end_vertex_id);
1063        if !at_corner {
1064            continue;
1065        }
1066        if boundary.contains(&edge.id) {
1067            if !(stations.contains_key(&edge.start_vertex_id)
1068                && stations.contains_key(&edge.end_vertex_id))
1069            {
1070                return Ok(None);
1071            }
1072            continue;
1073        }
1074        let incident = faces_of_edge.get(&edge.id).cloned().unwrap_or_default();
1075        if edge.degenerate
1076            || edge.curve.degree != 1
1077            || !incident.iter().all(|f| planar_neighbours.contains(f))
1078        {
1079            return Ok(None);
1080        }
1081    }
1082
1083    // The pushed radius, signed along the face's own outward normal.
1084    let [u0, u1] = face.surface.domain_u().or_refuse(KernelStage::Refine, "domain_u")?;
1085    let [v0, v1] = face.surface.domain_v().or_refuse(KernelStage::Refine, "domain_v")?;
1086    let (u_mid, v_mid) = (0.5 * (u0 + u1), 0.5 * (v0 + v1));
1087    let offsets = OffsetEvaluator::new(
1088        "push_revolution",
1089        &face.surface,
1090        OffsetNormal::Face {
1091            same_sense: face.same_sense,
1092        },
1093    );
1094    let sample = offsets.at(u_mid, v_mid, distance).or_refuse(KernelStage::Refine, "offset_at")?;
1095    let outward = cylindrical(sample.source).1.normalized().or_refuse(KernelStage::Refine, "normalized")?;
1096    let new_radius = cylindrical(sample.point).1.dot(outward);
1097    if new_radius <= tolerance {
1098        return Err(KernelRefusal::input(
1099            KernelStage::Refine,
1100            "radius_collapse",
1101            format!(
1102            "{op}: pushing partial cylinder face {face_id} by {distance} collapses its radius \
1103             {radius:.6} onto the axis — refusing"
1104        ),
1105        ));
1106    }
1107
1108    // Which meridian starts the sweep: the counter-clockwise run from start to
1109    // end is the one that passes through a rim's midpoint.
1110    let ccw = |from: f64, to: f64| (to - from).rem_euclid(TAU);
1111    let rim_mid = rims[0].0.curve.evaluate(0.5 * (rims[0].0.t0 + rims[0].0.t1)).or_refuse(KernelStage::Refine, "evaluate")?;
1112    let rim_azimuth = azimuth_of(cylindrical(rim_mid).1);
1113    let (start, end) = if ccw(meridians[0].azimuth, rim_azimuth)
1114        < ccw(meridians[0].azimuth, meridians[1].azimuth)
1115    {
1116        (&meridians[0], &meridians[1])
1117    } else {
1118        (&meridians[1], &meridians[0])
1119    };
1120
1121    // Solve one corner; `sense` is +1 where the arc leaves its corner
1122    // counter-clockwise (the start) and -1 where it leaves clockwise (the end).
1123    // Returns the new azimuth and the plane's two crossings of the new circle.
1124    let solve = |meridian: &Meridian, sense: f64| -> Result<(f64, [f64; 2]), KernelRefusal> {
1125        let normal = meridian.plane_normal;
1126        let (nx, ny) = (normal.dot(x_axis), normal.dot(y_axis));
1127        let reach = nx.hypot(ny);
1128        let phi = ny.atan2(nx);
1129        let offset = meridian.plane_origin.sub(origin).dot(normal);
1130        let lean = (phi - meridian.azimuth).sin();
1131        let side = if radius * reach * lean.abs() > residual_tolerance {
1132            sense * lean
1133        } else {
1134            -offset
1135        };
1136        let shortfall = offset.abs() - new_radius * reach;
1137        if shortfall > tolerance {
1138            return Err(KernelRefusal::unsupported(
1139                KernelStage::Refine,
1140                "partial_separated",
1141                format!(
1142                "{op}: pushing partial cylinder face {face_id} by {distance} separates it from \
1143                 neighbour face {} — the radius {new_radius:.6} carrier stops {shortfall:.3e} \
1144                 short of that plane (refusing)",
1145                meridian.neighbour
1146            ),
1147            ));
1148        }
1149        let half = (offset / (new_radius * reach)).clamp(-1.0, 1.0).acos();
1150        let azimuth = if sense * side > 0.0 { phi - half } else { phi + half };
1151        Ok((azimuth, [phi - half, phi + half]))
1152    };
1153    let (start_azimuth, start_crossings) = solve(start, 1.0)?;
1154    let (end_azimuth, end_crossings) = solve(end, -1.0)?;
1155    let new_sweep = ccw(start_azimuth, end_azimuth);
1156    let angular_tolerance = tolerance / new_radius;
1157    if new_sweep <= angular_tolerance || new_sweep >= TAU - angular_tolerance {
1158        return Err(KernelRefusal::input(
1159            KernelStage::Refine,
1160            "arc_closed",
1161            format!(
1162            "{op}: pushing partial cylinder face {face_id} by {distance} closes its arc between \
1163             neighbour faces {} and {} — refusing",
1164            start.neighbour, end.neighbour
1165        ),
1166        ));
1167    }
1168    // A crossing strictly inside the new arc means the arc passes through a
1169    // neighbour's plane before reaching its own corner — a convex round grown
1170    // past the two walls it was tangent to has no corner on either. Both
1171    // corners of an arc cut from one sketch line solve against two faces whose
1172    // planes agree only to the sketch solver's precision, so a crossing within
1173    // the residual bar of a corner is that corner.
1174    let crossing_margin = residual_tolerance / new_radius;
1175    for (meridian, crossings) in [(start, start_crossings), (end, end_crossings)] {
1176        for crossing in crossings {
1177            let along = ccw(start_azimuth, crossing);
1178            if along > crossing_margin && along < new_sweep - crossing_margin {
1179                return Err(KernelRefusal::unsupported(
1180                    KernelStage::Refine,
1181                    "arc_through_neighbour",
1182                    format!(
1183                    "{op}: pushing partial cylinder face {face_id} by {distance} carries its arc \
1184                     through the plane of neighbour face {} — refusing",
1185                    meridian.neighbour
1186                ),
1187                ));
1188            }
1189        }
1190    }
1191
1192    // New corner vertices, meridians and rims.
1193    let direction = |azimuth: f64| x_axis.scale(azimuth.cos()).add(y_axis.scale(azimuth.sin()));
1194    let mut new_vertices: HashMap<u64, Vec3> = HashMap::default();
1195    let mut new_curves: HashMap<u64, NurbsCurve> = HashMap::default();
1196    for (meridian, azimuth) in [(start, start_azimuth), (end, end_azimuth)] {
1197        let edge = meridian.edge;
1198        let mut ends = [Vec3::default(); 2];
1199        for (slot, vertex_id) in [edge.start_vertex_id, edge.end_vertex_id].into_iter().enumerate()
1200        {
1201            ends[slot] = origin
1202                .add(axis.scale(stations[&vertex_id]))
1203                .add(direction(azimuth).scale(new_radius));
1204            new_vertices.insert(vertex_id, ends[slot]);
1205        }
1206        new_curves.insert(edge.id, make_line(ends[0], ends[1]).or_refuse(KernelStage::Fragment, "make_line")?);
1207    }
1208    let new_x = direction(start_azimuth);
1209    let new_y = axis.cross(new_x);
1210    for (edge, station) in &rims {
1211        let (head, tail) = (new_vertices[&edge.start_vertex_id], new_vertices[&edge.end_vertex_id]);
1212        let mut arc = crate::make_arc(
1213            origin.add(axis.scale(*station)),
1214            new_x,
1215            new_y,
1216            new_radius,
1217            0.0,
1218            new_sweep,
1219        ).or_refuse(KernelStage::Fragment, "make_arc")?;
1220        let [a0, _] = arc.domain().or_refuse(KernelStage::Fragment, "domain")?;
1221        let first = arc.evaluate(a0).or_refuse(KernelStage::Fragment, "evaluate")?;
1222        if first.sub(head).length() > first.sub(tail).length() {
1223            arc = arc.reversed().or_refuse(KernelStage::Fragment, "reversed")?;
1224        }
1225        new_curves.insert(edge.id, arc);
1226    }
1227
1228    // The straight neighbour edges ending at a moved corner must still run the
1229    // way they ran; a corner carried past an edge's far end folds that face.
1230    for edge in &solid.edges {
1231        if new_curves.contains_key(&edge.id)
1232            || !(new_vertices.contains_key(&edge.start_vertex_id)
1233                || new_vertices.contains_key(&edge.end_vertex_id))
1234        {
1235            continue;
1236        }
1237        let position = |vertex_id: u64| {
1238            solid
1239                .vertices
1240                .iter()
1241                .find(|vertex| vertex.id == vertex_id)
1242                .map(|vertex| vertex.point)
1243                .ok_or_else(|| {
1244                    KernelRefusal::internal(
1245                        KernelStage::Fragment,
1246                        "corner_vertex_pos",
1247                        format!("{op}: missing vertex {vertex_id}"),
1248                    )
1249                })
1250        };
1251        let old_run = position(edge.end_vertex_id)?.sub(position(edge.start_vertex_id)?);
1252        let new_start = match new_vertices.get(&edge.start_vertex_id) {
1253            Some(point) => *point,
1254            None => position(edge.start_vertex_id)?,
1255        };
1256        let new_end = match new_vertices.get(&edge.end_vertex_id) {
1257            Some(point) => *point,
1258            None => position(edge.end_vertex_id)?,
1259        };
1260        let new_run = new_end.sub(new_start);
1261        if new_run.length() <= tolerance || new_run.dot(old_run) <= 0.0 {
1262            return Err(KernelRefusal::unsupported(
1263                KernelStage::Fragment,
1264                "corner_past_edge",
1265                format!(
1266                "{op}: pushing partial cylinder face {face_id} by {distance} runs a corner past \
1267                 the far end of edge {} — refusing",
1268                edge.id
1269            ),
1270            ));
1271        }
1272    }
1273
1274    // The offset carrier: the same axis, the new radius, the new sweep, with
1275    // the generatrix kept in its own direction so the normal keeps its sense.
1276    let carrier = make_revolution(
1277        origin,
1278        axis,
1279        &make_line(
1280            origin
1281                .add(axis.scale(generatrix_z0))
1282                .add(new_x.scale(new_radius)),
1283            origin
1284                .add(axis.scale(generatrix_z1))
1285                .add(new_x.scale(new_radius)),
1286        ).or_refuse(KernelStage::Fragment, "make_line")?,
1287        new_sweep,
1288    ).or_refuse(KernelStage::Fragment, "make_revolution")?;
1289    let [c0, c1] = carrier.domain_u().or_refuse(KernelStage::Fragment, "domain_u")?;
1290    let [w0, w1] = carrier.domain_v().or_refuse(KernelStage::Fragment, "domain_v")?;
1291    let old_normal = face.surface.normal(u_mid, v_mid).or_refuse(KernelStage::Fragment, "normal")?;
1292    let old_point = face.surface.evaluate(u_mid, v_mid).or_refuse(KernelStage::Fragment, "evaluate")?;
1293    let new_normal = carrier.normal(0.5 * (c0 + c1), 0.5 * (w0 + w1)).or_refuse(KernelStage::Fragment, "normal")?;
1294    let new_point = carrier.evaluate(0.5 * (c0 + c1), 0.5 * (w0 + w1)).or_refuse(KernelStage::Fragment, "evaluate")?;
1295    if old_normal.dot(cylindrical(old_point).1) * new_normal.dot(cylindrical(new_point).1) <= 0.0 {
1296        return Err(KernelRefusal::internal(
1297            KernelStage::Refine,
1298            "carrier_normal",
1299            format!(
1300            "{op}: the offset carrier of partial cylinder face {face_id} flipped its normal — \
1301             refusing"
1302        ),
1303        ));
1304    }
1305
1306    let mut result = solid.clone();
1307    let pushed = &mut result.shells[shell_index].faces[face_index];
1308    for coedge in &mut pushed.loops[0].coedges {
1309        let mut pcurve = build_pcurve_on_surface(&carrier, &new_curves[&coedge.edge_id]).or_refuse(KernelStage::Sew, "build_pcurve_on_surface")?;
1310        if !coedge.forward {
1311            pcurve = pcurve.reversed().or_refuse(KernelStage::Sew, "reversed")?;
1312        }
1313        coedge.pcurve = pcurve;
1314    }
1315    pushed.surface = carrier;
1316    retrim_planar_neighbours(
1317        solid,
1318        result,
1319        &edge_by_id,
1320        &new_curves,
1321        &new_vertices,
1322        planar_neighbours,
1323        residual_tolerance,
1324        scale,
1325    )
1326    .map(Some)
1327}
1328
1329/// Land a push's rebuilt boundary on `result` (whose pushed face already
1330/// carries its new surface and pcurves): swap in the new edge curves and vertex
1331/// positions, rebuild every straight edge of a planar neighbour that ends at a
1332/// moved vertex, re-trim those neighbours, and accept the solid only if it
1333/// validates and keeps its orientation.
1334#[allow(clippy::too_many_arguments)]
1335fn retrim_planar_neighbours(
1336    solid: &BrepSolid,
1337    mut result: BrepSolid,
1338    edge_by_id: &HashMap<u64, &EdgeRecord>,
1339    new_curves: &HashMap<u64, NurbsCurve>,
1340    new_vertices: &HashMap<u64, Vec3>,
1341    planar_neighbours: HashSet<u64>,
1342    residual_tolerance: f64,
1343    scale: f64,
1344) -> Result<BrepSolid, KernelRefusal> {
1345    for edge in &mut result.edges {
1346        if let Some(curve) = new_curves.get(&edge.id) {
1347            let [d0, d1] = curve.domain().or_refuse(KernelStage::Sew, "domain")?;
1348            edge.curve = curve.clone();
1349            edge.t0 = d0;
1350            edge.t1 = d1;
1351        }
1352    }
1353    for vertex in &mut result.vertices {
1354        if let Some(point) = new_vertices.get(&vertex.id) {
1355            vertex.point = *point;
1356        }
1357    }
1358    // A full-revolve planar disk/annulus may carry its own radial seam from an
1359    // axis vertex to the shared circular rim. Moving the rim relocates one end
1360    // of that seam, so rebuild the straight edge before re-trimming the cap.
1361    let result_vertices: HashMap<u64, Vec3> = result
1362        .vertices
1363        .iter()
1364        .map(|vertex| (vertex.id, vertex.point))
1365        .collect();
1366    let mut cap_seams = HashSet::default();
1367    for neighbour in &planar_neighbours {
1368        let (ns, nf) = find_face(&result, *neighbour).ok_or_else(|| {
1369            KernelRefusal::internal(
1370                KernelStage::Sew,
1371                "cap_lookup",
1372                format!("offset_revolution_face: missing cap {neighbour}"),
1373            )
1374        })?;
1375        for loop_record in &result.shells[ns].faces[nf].loops {
1376            for coedge in &loop_record.coedges {
1377                let edge = *edge_by_id.get(&coedge.edge_id).ok_or_else(|| {
1378                    KernelRefusal::internal(
1379                        KernelStage::Sew,
1380                        "missing_cap_edge",
1381                        format!(
1382                        "offset_revolution_face: missing cap edge {}",
1383                        coedge.edge_id
1384                    ),
1385                    )
1386                })?;
1387                let touches_moved = new_vertices.contains_key(&edge.start_vertex_id)
1388                    || new_vertices.contains_key(&edge.end_vertex_id);
1389                if touches_moved
1390                    && !new_curves.contains_key(&edge.id)
1391                    && !edge.degenerate
1392                    && edge.curve.degree == 1
1393                {
1394                    cap_seams.insert(edge.id);
1395                }
1396            }
1397        }
1398    }
1399    for seam_id in cap_seams {
1400        let source = edge_by_id[&seam_id];
1401        let curve = make_line(
1402            result_vertices[&source.start_vertex_id],
1403            result_vertices[&source.end_vertex_id],
1404        ).or_refuse(KernelStage::Sew, "make_line")?;
1405        let [d0, d1] = curve.domain().or_refuse(KernelStage::Sew, "domain")?;
1406        if let Some(edge) = result.edges.iter_mut().find(|edge| edge.id == seam_id) {
1407            edge.curve = curve;
1408            edge.t0 = d0;
1409            edge.t1 = d1;
1410        }
1411    }
1412    let final_edges: HashMap<u64, EdgeRecord> = result
1413        .edges
1414        .iter()
1415        .map(|edge| (edge.id, edge.clone()))
1416        .collect();
1417    for neighbour in planar_neighbours {
1418        let (ns, nf) = find_face(&result, neighbour).ok_or_else(|| {
1419            KernelRefusal::internal(
1420                KernelStage::Sew,
1421                "retrim_cap",
1422                format!("offset_revolution_face: missing cap {neighbour}"),
1423            )
1424        })?;
1425        let plane = plane_of_surface(
1426            &result.shells[ns].faces[nf].surface,
1427            residual_tolerance,
1428            "offset_revolution_face",
1429        )?;
1430        retrim_planar_face(
1431            &mut result.shells[ns].faces[nf],
1432            &plane,
1433            &final_edges,
1434            scale,
1435            "offset_revolution_face",
1436        )
1437        .map_err(|error| {
1438            KernelRefusal::internal(
1439                KernelStage::Sew,
1440                "cap_retrim",
1441                format!("offset_revolution_face: cap retrim: {error}"),
1442            )
1443        })?;
1444    }
1445    let issues = result.validate();
1446    if !issues.is_empty() {
1447        return Err(KernelRefusal::internal(
1448            KernelStage::Validate,
1449            "validate",
1450            format!(
1451            "offset_revolution_face: pushed solid failed validation: {issues:?}"
1452        ),
1453        ));
1454    }
1455    if let (Ok(before), Ok(after)) = (solid_signed_volume(solid), solid_signed_volume(&result)) {
1456        if before * after <= 0.0 {
1457            return Err(KernelRefusal::internal(
1458                KernelStage::Validate,
1459                INVERTED,
1460                "offset_revolution_face: the push inverts the solid — refusing",
1461            ));
1462        }
1463    }
1464    Ok(result)
1465}
1466
1467/// `BREP_OFFSET_REVOLUTION_TRACE=<file>` appends one record per general-lane
1468/// push to `<file>`: the call, each rebuilt boundary edge against the fixed
1469/// planar neighbour it must lie in — before the push (the source edge) and
1470/// after (the rebuilt curve) — and the outcome with the volume. It writes to a
1471/// file, not stderr, because the case gate replays each document in a child
1472/// process whose stderr it discards.
1473struct BoundaryTrace {
1474    path: Option<std::path::PathBuf>,
1475}
1476
1477impl BoundaryTrace {
1478    fn open(solid: &BrepSolid, face_id: u64, distance: f64) -> Self {
1479        let path = std::env::var_os("BREP_OFFSET_REVOLUTION_TRACE")
1480            .filter(|value| !value.is_empty() && value != "0")
1481            .map(std::path::PathBuf::from);
1482        let trace = Self { path };
1483        if trace.path.is_some() {
1484            let volume = solid_signed_volume(solid).map_or(f64::NAN, f64::abs);
1485            trace.write(format!(
1486                "call\tface={face_id}\tdistance={distance}\tscale={}\tvolume_before={volume:.15}",
1487                solid_model_scale(solid)
1488            ));
1489        }
1490        trace
1491    }
1492
1493    /// One whole line per `write` call on an append-mode file, prefixed with
1494    /// the thread's name, so parallel test threads (named after their tests)
1495    /// neither interleave nor lose which test a line belongs to.
1496    fn write(&self, line: String) {
1497        use std::io::Write;
1498        let Some(path) = &self.path else { return };
1499        let thread = std::thread::current();
1500        let record = format!("{}\t{line}\n", thread.name().unwrap_or("?"));
1501        if let Ok(mut file) = std::fs::OpenOptions::new().create(true).append(true).open(path) {
1502            let _ = file.write_all(record.as_bytes());
1503        }
1504    }
1505
1506    /// Largest distance from `plane` over `samples + 1` points of `curve`
1507    /// between `t0` and `t1`, and over its two end points.
1508    fn off_plane(
1509        curve: &NurbsCurve,
1510        t0: f64,
1511        t1: f64,
1512        plane: &crate::offset_retrim::Plane,
1513    ) -> (f64, f64) {
1514        let distance = |t: f64| {
1515            curve
1516                .evaluate(t)
1517                .map_or(f64::NAN, |point| point.sub(plane.origin).dot(plane.normal).abs())
1518        };
1519        let ends = distance(t0).max(distance(t1));
1520        let samples = (0..=32)
1521            .map(|index| distance(t0 + (t1 - t0) * index as f64 / 32.0))
1522            .fold(0.0f64, f64::max);
1523        (ends, samples)
1524    }
1525
1526    fn edge(
1527        &self,
1528        kind: &str,
1529        edge: &EdgeRecord,
1530        neighbour: Option<&u64>,
1531        plane: Option<&crate::offset_retrim::Plane>,
1532        rebuilt: &NurbsCurve,
1533    ) {
1534        if self.path.is_none() {
1535            return;
1536        }
1537        let Some(plane) = plane else {
1538            self.write(format!("edge\t{}\t{kind}\tneighbour=none", edge.id));
1539            return;
1540        };
1541        let (before_ends, before_samples) = Self::off_plane(&edge.curve, edge.t0, edge.t1, plane);
1542        let (after_ends, after_samples) = match rebuilt.domain() {
1543            Ok([r0, r1]) => Self::off_plane(rebuilt, r0, r1, plane),
1544            Err(_) => (f64::NAN, f64::NAN),
1545        };
1546        self.write(format!(
1547            "edge\t{}\t{kind}\tneighbour={}\tbefore_vertices={before_ends:.3e}\t\
1548             before_samples={before_samples:.3e}\tafter_vertices={after_ends:.3e}\t\
1549             after_samples={after_samples:.3e}",
1550            edge.id,
1551            neighbour.map_or("none".to_string(), u64::to_string)
1552        ));
1553    }
1554
1555    fn outcome(&self, outcome: &Result<BrepSolid, KernelRefusal>) {
1556        if self.path.is_none() {
1557            return;
1558        }
1559        match outcome {
1560            Ok(result) => {
1561                let volume = solid_signed_volume(result).map_or(f64::NAN, f64::abs);
1562                self.write(format!("ok\tvolume_after={volume:.15}"));
1563            }
1564            Err(error) => self.write(format!("refused\t{error}")),
1565        }
1566    }
1567}
1568