Skip to main content

ogeom_offset/
draft.rs

1//! Draft: turning faces about a neutral plane so a part can leave its mould.
2//!
3//! A drafted face is the same face on a *tilted* support. It keeps the line
4//! where it crosses the neutral plane (that line does not move, which is
5//! what makes the draft measurable from a datum) and turns about it by the
6//! draft angle. Everything else follows: the neighbouring faces re-meet the
7//! tilted plane, the vertices re-solve, and the solid comes back with the
8//! same topology on new geometry.
9//!
10//! That last part is not this module's work. It is the offset's rebuild,
11//! which puts a solid back together on moved supports. A draft hands it
12//! turned surfaces instead of translated ones.
13
14use ogeom_algo::Built;
15use ogeom_core::{OgeomResult, Tolerances, ogeom_bail};
16use ogeom_geom::{PlaneSurface, Surface as _, SurfaceGeometry};
17use ogeom_math::{Direction, Frame, Plane, Point, Transform, Vector};
18use ogeom_topo::{Model, NodeData, Shape, ShapeType, TShapeId};
19
20use crate::shape::rebuilt;
21
22/// Draft the named faces of a solid about a neutral plane.
23///
24/// Each face turns about its own intersection with `neutral` by `angle`,
25/// in the sense that leans the face inwards as it goes: a positive angle
26/// narrows the solid in the `pull` direction (the way the part leaves its
27/// mould), and a negative one widens it. Leaning inwards tilts the face's
28/// outward normal *towards* the pull, which is how the sense is picked, by
29/// measurement. A face the neutral plane does not reach turns about where
30/// its surface, run on beyond the face, meets the plane. A face parallel
31/// to the neutral plane has no line to turn about and is refused by name,
32/// as is a face the plane crosses in two runs (a region beside both would
33/// turn about two hinges), a patch that stops short of the plane, and a
34/// face the rebuild cannot re-meet.
35///
36/// # Errors
37///
38/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if a
39/// named face is not a planar face of `solid`, is parallel to the neutral
40/// plane, or the angle is not a usable one; plus whatever the rebuild
41/// refuses.
42pub fn apply_draft(
43    model: &mut Model,
44    solid: &Shape,
45    faces: &[Shape],
46    neutral: Plane,
47    pull: Direction,
48    angle: f64,
49    tol: Tolerances,
50) -> OgeomResult<Built> {
51    if !angle.is_finite() || angle.abs() >= core::f64::consts::FRAC_PI_2 {
52        ogeom_bail!(
53            Construction,
54            "a draft of {angle} radians turns the face past its own plane"
55        );
56    }
57    let (canonical, mapped, prefix) = crate::shape::canonical_input(model, solid, faces, tol)?;
58    if let Some(prefix) = prefix {
59        let mut out = apply_draft(model, &canonical, &mapped, neutral, pull, angle, tol)?;
60        out.history = prefix.then(&out.history);
61        return Ok(out);
62    }
63    if faces.is_empty() {
64        ogeom_bail!(Construction, "a draft of no faces drafts nothing");
65    }
66    // The solid's own face occurrences, orientation and all: a handle a
67    // caller got from a canonical exploration carries no use-orientation,
68    // and the draft's sense probe needs the true outward.
69    let own: Vec<Shape> = {
70        let mut seen: Vec<Shape> = Vec::new();
71        for f in ogeom_topo::explore(model, solid, ogeom_topo::Filter::OfType(ShapeType::Face))? {
72            if !seen.iter().any(|s| s.node() == f.node()) {
73                seen.push(f);
74            }
75        }
76        seen
77    };
78
79    // The turned surface for each named face, worked out before the
80    // rebuild, so a face that cannot be drafted says so here rather than
81    // half-way through a solid.
82    let mut turned: Vec<(TShapeId, SurfaceGeometry)> = Vec::with_capacity(faces.len());
83    for face in faces {
84        let Some(used) = own.iter().find(|f| f.node() == face.node()).cloned() else {
85            ogeom_bail!(Construction, "a drafted face is not a face of the solid");
86        };
87        let face = &used;
88        let Some(NodeData::Face(data)) = model.node(face).map(|n| n.data().clone()) else {
89            ogeom_bail!(Construction, "expected a face");
90        };
91        let Some(surface) = model.geometry().surface(data.surface) else {
92            ogeom_bail!(Dangling, "face refers to a surface not in this model");
93        };
94        // Which sign turns the raw surface normal outward, read from the
95        // material itself rather than from an orientation flag a handle may
96        // or may not carry: a step along the raw normal that lands inside
97        // the solid means the raw normal points inward.
98        let sign = outward_sign(model, solid, face, surface, tol)?;
99        let sign_of = |_: &Shape| sign;
100        // A wall of revolution drafts about its neutral *circle*: the same
101        // axis, the radius at the neutral plane held, the slant turned.
102        let axial = |frame: Frame| -> bool {
103            (frame.z().vector().dot(neutral.normal().vector()).abs() - 1.0).abs()
104                <= tol.angular().max(1e-9)
105        };
106        match surface {
107            SurfaceGeometry::Cylinder(c) if axial(c.cylinder().frame()) => {
108                let cylinder = c.cylinder();
109                let (_, (v0, v1)) = surface.domain();
110                turned.push((
111                    face.node(),
112                    revolved_draft(
113                        cylinder.frame(),
114                        cylinder.radius(),
115                        0.0,
116                        (v0, v1),
117                        sign_of(face),
118                        neutral,
119                        pull,
120                        angle,
121                        tol,
122                    )?,
123                ));
124                continue;
125            }
126            SurfaceGeometry::Cone(co) if axial(co.cone().frame()) => {
127                let cone = co.cone();
128                let (_, (v0, v1)) = surface.domain();
129                turned.push((
130                    face.node(),
131                    revolved_draft(
132                        cone.frame(),
133                        cone.reference_radius(),
134                        cone.half_angle(),
135                        (v0, v1),
136                        sign_of(face),
137                        neutral,
138                        pull,
139                        angle,
140                        tol,
141                    )?,
142                ));
143                continue;
144            }
145            SurfaceGeometry::Extrusion(e) => {
146                turned.push((
147                    face.node(),
148                    extruded_draft(
149                        e,
150                        surface.domain(),
151                        sign_of(face),
152                        neutral,
153                        pull,
154                        angle,
155                        tol,
156                    )?,
157                ));
158                continue;
159            }
160            SurfaceGeometry::Plane(_) => {}
161            // Everything else (a raw fitted patch, a wall of revolution
162            // about an oblique neutral) is drafted the way a mould-maker
163            // drafts: along the pull, tilted by the angle, through the
164            // line where the face crosses the neutral plane.
165            _ => {
166                turned.push((
167                    face.node(),
168                    general_draft(model, face, surface, sign, neutral, pull, angle, tol)?,
169                ));
170                continue;
171            }
172        }
173        let SurfaceGeometry::Plane(p) = surface else {
174            unreachable!("the match above let only planes through");
175        };
176        let plane = p.plane();
177        let ((u0, u1), (v0, v1)) = surface.domain();
178        // The plane's raw normal is its du x dv; the measured sign turns it
179        // outward.
180        let outward = plane.normal().vector() * sign;
181
182        // The hinge: the line where this face crosses the neutral plane.
183        let along = plane.normal().vector().cross(neutral.normal().vector());
184        let magnitude = along.magnitude();
185        if magnitude <= tol.angular() {
186            ogeom_bail!(
187                Construction,
188                "a face parallel to the neutral plane has no line to turn \
189                 about"
190            );
191        }
192        let along = along / magnitude;
193        let hinge = meet(plane, neutral, along, tol)?;
194
195        // Which way to turn: probed at the angle's *magnitude*, so the
196        // sense names the inward lean (outward normal furthest towards the
197        // pull, the solid narrowing as it leaves), and the angle's sign
198        // stays the caller's: positive drafts inward, negative outward.
199        let axis = ogeom_math::Axis::new(hinge, Direction::new(along, tol)?);
200        let mut candidates = Vec::with_capacity(2);
201        for sense in [1.0, -1.0] {
202            let turn = Transform::rotation(axis, angle.abs() * sense);
203            candidates.push((sense, turn.apply_vector(outward).dot(pull.vector())));
204        }
205        let leaning = candidates
206            .iter()
207            .copied()
208            .max_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(core::cmp::Ordering::Equal))
209            .map_or(1.0, |(sense, _)| sense);
210        let turn = Transform::rotation(axis, angle * leaning);
211        let moved_normal = Direction::new(turn.apply_vector(plane.normal().vector()), tol)?;
212        let tilted = Plane::new(Frame::new(
213            hinge,
214            moved_normal,
215            Direction::new(along, tol)?,
216            tol,
217        )?);
218        // The old window turned with the face, read in the tilted plane's
219        // own chart: its origin is the hinge and its `u` runs along it, so
220        // the old window's numbers name another stretch of it.
221        let ((u0, u1), (v0, v1)) = if [u0, u1, v0, v1].iter().all(|x| x.is_finite()) {
222            let (mut lo, mut hi) = (
223                (f64::INFINITY, f64::INFINITY),
224                (f64::NEG_INFINITY, f64::NEG_INFINITY),
225            );
226            for (u, v) in [(u0, v0), (u1, v0), (u1, v1), (u0, v1)] {
227                let corner = tilted
228                    .frame()
229                    .to_local(turn.apply(surface.point_at(u, v, tol)?));
230                lo = (lo.0.min(corner.x), lo.1.min(corner.y));
231                hi = (hi.0.max(corner.x), hi.1.max(corner.y));
232            }
233            ((lo.0, hi.0), (lo.1, hi.1))
234        } else {
235            ((u0, u1), (v0, v1))
236        };
237        // The window grows with the turn: a tilted plane reaches further
238        // across the same solid than the one it replaces.
239        let grow = (u1 - u0).abs().max((v1 - v0).abs()).mul_add(0.5, 1.0) * angle.abs().tan()
240            + tol.confusion();
241        turned.push((
242            face.node(),
243            PlaneSurface::over(tilted, (u0 - grow, u1 + grow), (v0 - grow, v1 + grow))?.into(),
244        ));
245    }
246
247    rebuilt(
248        model,
249        solid,
250        &|_| 0.0,
251        &|face| {
252            turned
253                .iter()
254                .find(|(node, _)| *node == face.node())
255                .map(|(_, surface)| surface.clone())
256        },
257        tol,
258    )
259}
260
261/// The turned support for a drafted wall of revolution: a cone about the
262/// same axis, holding the radius at the neutral plane and leaning the slant
263/// by the draft, in the sense that tips the outward normal towards the pull.
264#[allow(clippy::too_many_arguments, reason = "one construction, all its data")]
265fn revolved_draft(
266    frame: Frame,
267    reference_radius: f64,
268    half_angle: f64,
269    window: (f64, f64),
270    sign: f64,
271    neutral: Plane,
272    pull: Direction,
273    angle: f64,
274    tol: Tolerances,
275) -> OgeomResult<SurfaceGeometry> {
276    use ogeom_geom::ConeSurface;
277
278    let axis_dir = frame.z().vector();
279    let along = axis_dir.dot(neutral.normal().vector());
280    if (along.abs() - 1.0).abs() > tol.angular().max(1e-9) {
281        ogeom_bail!(
282            Construction,
283            "a wall of revolution drafts about a neutral plane square to \
284             its axis; the oblique neutral needs the general machinery (see \
285             docs/CAPABILITIES.md, offset.draft)"
286        );
287    }
288    // The neutral circle: where the axis meets the plane, and the radius
289    // the wall holds there.
290    let height = -neutral.signed_distance_to(frame.origin()) * along.signum();
291    let neutral_point = frame.origin() + axis_dir * height;
292    let neutral_radius = half_angle.tan().mul_add(height, reference_radius);
293    if neutral_radius <= tol.confusion() {
294        ogeom_bail!(
295            Construction,
296            "the wall has no radius left at the neutral plane to hold"
297        );
298    }
299    let hinge_frame = Frame::new(neutral_point, frame.z(), frame.x(), tol)?;
300
301    // Which way to lean, by measurement: of the two candidate slants, keep
302    // the one whose outward normal (probed a little above the neutral
303    // circle) ends up leaning furthest towards the pull.
304    let mut best: Option<(f64, f64)> = None;
305    for sense in [1.0_f64, -1.0] {
306        // Probed at the magnitude: the sense names the inward lean, and the
307        // caller's sign then picks inward or outward through it.
308        let probe = half_angle + angle.abs() * sense;
309        let candidate = half_angle + angle * sense;
310        if probe.abs() <= tol.angular()
311            || probe.abs() >= core::f64::consts::FRAC_PI_2 - tol.angular()
312            || candidate.abs() <= tol.angular()
313            || candidate.abs() >= core::f64::consts::FRAC_PI_2 - tol.angular()
314        {
315            continue;
316        }
317        let cone = ogeom_math::Cone::new(hinge_frame, neutral_radius, probe, tol)?;
318        let surface: SurfaceGeometry = ConeSurface::new(cone, (-1.0, 1.0))?.into();
319        let (du, dv) = surface.d1_at(0.0, 1.0, tol)?;
320        let n = du.cross(dv);
321        let outward = n / n.magnitude() * sign;
322        let lean = outward.dot(pull.vector());
323        if best.as_ref().is_none_or(|(_, held)| lean > *held) {
324            best = Some((candidate, lean));
325        }
326    }
327    let Some((leaned, _)) = best else {
328        ogeom_bail!(
329            Construction,
330            "a draft of {angle} radians flattens the wall or swallows it"
331        );
332    };
333    let cone = ogeom_math::Cone::new(hinge_frame, neutral_radius, leaned, tol)?;
334
335    // The old window, re-expressed against the neutral origin and grown a
336    // little; refused when the slant runs out of radius inside it.
337    let shift = height;
338    let grow = (window.1 - window.0).abs().mul_add(0.1, 1.0);
339    let (w0, w1) = (window.0 - shift - grow, window.1 - shift + grow);
340    let apex_height = -neutral_radius / leaned.tan();
341    if apex_height > w0 && apex_height < w1 {
342        ogeom_bail!(
343            Construction,
344            "the draft swallows the drafted face's own apex"
345        );
346    }
347    Ok(ConeSurface::new(cone, (w0, w1))?.into())
348}
349
350/// The turned support for a drafted extruded wall: every ruling rotated
351/// about the hinge curve's own tangent by the draft, the result re-fitted.
352///
353/// The hinge is where the wall crosses the neutral plane (one closed-form
354/// height per profile parameter), and it does not move, exactly as a planar
355/// draft's hinge line does not. Each ruling turns about the hinge's local
356/// tangent in the sense that leans the outward normal towards the pull,
357/// probed at the profile's midpoint the way the planar draft probes its
358/// candidates. The turned rulings are sampled on a grid and fitted; a draft
359/// whose rulings cross inside the drafted window (a concave profile turned
360/// far enough to fold) is refused by name before anything is fitted.
361#[allow(clippy::too_many_arguments, reason = "one construction, all its data")]
362fn extruded_draft(
363    extrusion: &ogeom_geom::ExtrusionSurface,
364    window: ((f64, f64), (f64, f64)),
365    sign: f64,
366    neutral: Plane,
367    pull: Direction,
368    angle: f64,
369    tol: Tolerances,
370) -> OgeomResult<SurfaceGeometry> {
371    use ogeom_geom::Curve3d as _;
372    let ((u0, u1), (v0, v1)) = window;
373    let d = extrusion.direction().vector();
374    let n = neutral.normal().vector();
375    let den = n.dot(d);
376    if den.abs() <= tol.angular() {
377        ogeom_bail!(
378            Construction,
379            "the neutral plane runs along the wall's rulings; there is no \
380             hinge to turn about"
381        );
382    }
383    let curve = extrusion.curve();
384    let o = neutral.origin().to_vector();
385    // Where the ruling through C(u) crosses the neutral plane, and which way
386    // the hinge runs there.
387    let height_at = |c: Point| n.dot(o - c.to_vector()) / den;
388    let hinge_tangent = |cd: Vector| cd - d * (n.dot(cd) / den);
389
390    // The sense, probed at the profile's midpoint exactly as the planar
391    // draft probes its two candidates: the turn whose outward normal leans
392    // furthest towards the pull is the inward one, and the caller's sign
393    // picks inward or outward through it.
394    let um = f64::midpoint(u0, u1);
395    let cm = curve.point_at(um, tol)?;
396    let cdm = curve.d1_at(um, tol)?;
397    let hinge_m = cm + d * height_at(cm);
398    let tangent_m = Direction::new(hinge_tangent(cdm), tol)?;
399    let outward = {
400        let nw = cdm.cross(d);
401        nw / nw.magnitude() * sign
402    };
403    let axis_m = ogeom_math::Axis::new(hinge_m, tangent_m);
404    let mut leaning = 1.0;
405    let mut best = f64::NEG_INFINITY;
406    for sense in [1.0_f64, -1.0] {
407        let turn = Transform::rotation(axis_m, angle.abs() * sense);
408        let lean = turn.apply_vector(outward).dot(pull.vector());
409        if lean > best {
410            best = lean;
411            leaning = sense;
412        }
413    }
414    let theta = angle * leaning;
415
416    // One ruling per profile parameter: the hinge point, the hinge
417    // tangent, and the extrusion direction turned about it.
418    let ruling_at = |u: f64| -> OgeomResult<(Point, Vector, f64)> {
419        let c = curve.point_at(u, tol)?;
420        let cd = curve.d1_at(u, tol)?;
421        let h = height_at(c);
422        let hinge = c + d * h;
423        let tangent = Direction::new(hinge_tangent(cd), tol)?;
424        let turn = Transform::rotation(ogeom_math::Axis::new(hinge, tangent), theta);
425        Ok((hinge, turn.apply_vector(d), h))
426    };
427    const ALONG: usize = 65;
428    #[allow(clippy::cast_precision_loss)]
429    let along_u = |x: f64| u0 + (u1 - u0) * x / ((ALONG - 1) as f64);
430    let mut base: Vec<(Point, Vector)> = Vec::with_capacity(ALONG);
431    let (mut s_lo, mut s_hi) = (f64::INFINITY, f64::NEG_INFINITY);
432    for i in 0..ALONG {
433        #[allow(clippy::cast_precision_loss)]
434        let (hinge, ruling, h) = ruling_at(along_u(i as f64))?;
435        base.push((hinge, ruling));
436        s_lo = s_lo.min(v0 - h);
437        s_hi = s_hi.max(v1 - h);
438    }
439    // The window grows with the turn, as the planar draft's does. Along the
440    // profile the curve itself ends, so the growth is a tangent-line
441    // continuation at each end: the wall must still reach the neighbours it
442    // re-meets, and the continuation exists only to be trimmed away there.
443    let grow = (u1 - u0).abs().max((v1 - v0).abs()).mul_add(0.5, 1.0) * angle.abs().tan()
444        + tol.confusion();
445    let (s_lo, s_hi) = (s_lo - grow, s_hi + grow);
446    // Quadratic continuation, so the fitted wall keeps its end curvature
447    // across the join instead of kinking straight: `k` station spacings
448    // past the end, from the end's last three stations.
449    struct Continuation {
450        hinge: Point,
451        d1: Vector,
452        d2: Vector,
453        ruling: Vector,
454        r1: Vector,
455        steps: usize,
456    }
457    let continuation = |front: bool| -> Continuation {
458        let (i0, i1, i2) = if front {
459            (0, 1, 2)
460        } else {
461            (ALONG - 1, ALONG - 2, ALONG - 3)
462        };
463        let d1 = base[i0].0 - base[i1].0;
464        let d2 = (base[i0].0 - base[i1].0) - (base[i1].0 - base[i2].0);
465        // Far enough past the end that a neighbour standing across the
466        // wall at a slant still finds wall where the turned rulings have
467        // carried it sideways, as far as they lean at the window's far
468        // heights.
469        let step = d1.magnitude().max(tol.confusion());
470        let ruling = base[i0].1;
471        let lean = (ruling - d * ruling.dot(d)).magnitude() * s_lo.abs().max(s_hi.abs());
472        let steps = ((grow + lean) / step).ceil().max(2.0);
473        #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
474        let steps = (steps as usize).min(16);
475        Continuation {
476            hinge: base[i0].0,
477            d1,
478            d2,
479            ruling: base[i0].1,
480            r1: base[i0].1 - base[i1].1,
481            steps,
482        }
483    };
484    let (front, back) = (continuation(true), continuation(false));
485    let continued = |c: &Continuation, k: f64| -> (Point, Vector) {
486        (
487            c.hinge + c.d1 * k + c.d2 * (k * (k + 1.0) / 2.0),
488            c.ruling + c.r1 * k,
489        )
490    };
491    // The wall's stations by index `x` along it: the front continuation,
492    // the profile's stations, the back continuation; between the
493    // profile's stations the profile itself, and between the
494    // continuation's its formula.
495    #[allow(clippy::cast_precision_loss)]
496    let (first, last) = (front.steps as f64, (front.steps + ALONG - 1) as f64);
497    let station = |x: f64| -> OgeomResult<(Point, Vector)> {
498        if x < first {
499            return Ok(continued(&front, first - x));
500        }
501        if x > last {
502            return Ok(continued(&back, x - last));
503        }
504        let (hinge, ruling, _) = ruling_at(along_u(x - first))?;
505        Ok((hinge, ruling))
506    };
507    let along_total = front.steps + ALONG + back.steps;
508    let stations: Vec<(Point, Vector)> = (0..along_total)
509        .map(|i| {
510            #[allow(clippy::cast_precision_loss)]
511            let x = i as f64;
512            match i.checked_sub(front.steps) {
513                Some(j) if j < ALONG => Ok(base[j]),
514                _ => station(x),
515            }
516        })
517        .collect::<OgeomResult<_>>()?;
518
519    // A fold is two rulings crossing inside the window: walking the wall at
520    // either extreme height must still advance the way the hinge advances.
521    for edge in [s_lo, s_hi] {
522        for pair in stations.windows(2) {
523            let ((h0, r0), (h1, r1)) = (pair[0], pair[1]);
524            let step = (h1 + r1 * edge) - (h0 + r0 * edge);
525            if step.dot(h1 - h0) <= 0.0 {
526                ogeom_bail!(
527                    Construction,
528                    "the draft folds the wall onto itself inside the drafted \
529                     window; refused; see docs/CAPABILITIES.md, offset.draft"
530                );
531            }
532        }
533    }
534
535    // Rulings are straight, so a handful of rows fits them exactly; the
536    // profile direction carries the shape, and is fitted to the turned
537    // rulings at the stations and between them, the stations refined
538    // where the wall misses.
539    const ACROSS: usize = 9;
540    // The chart runs over [0, 1] both ways: `u` across the window from its
541    // top down, `v` along the wall by station index, which keeps the
542    // extrusion's own sense (`du x dv` the same way round). The straight
543    // rulings take `u`, where a measure of the face integrates along it.
544    #[allow(clippy::cast_precision_loss)]
545    let span = (along_total - 1) as f64;
546    let xs: Vec<f64> = (0..along_total)
547        .map(|i| {
548            #[allow(clippy::cast_precision_loss)]
549            let x = i as f64 / span;
550            x
551        })
552        .collect();
553    let ss: Vec<f64> = (0..ACROSS)
554        .map(|j| {
555            #[allow(clippy::cast_precision_loss)]
556            let s = j as f64 / ((ACROSS - 1) as f64);
557            s
558        })
559        .collect();
560    let fit_target = (tol.confusion() * 1e3).max(1e-4);
561    let fitted = ogeom_geom::fit::fit_surface_sampled(
562        |s, x| {
563            let x = x * span;
564            #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
565            let i = x.round() as usize;
566            #[allow(clippy::cast_precision_loss)]
567            let (hinge, ruling) = if (i as f64 - x).abs() <= 1e-12 * span && i < along_total {
568                stations[i]
569            } else {
570                station(x)?
571            };
572            Ok(hinge + ruling * (s_hi - (s_hi - s_lo) * s))
573        },
574        &ogeom_geom::fit::Sampling {
575            us: ss,
576            vs: xs,
577            between: (false, true),
578            closed_v: false,
579            most: 4096,
580        },
581        3,
582        fit_target,
583        tol,
584    )?;
585    if !fitted.met {
586        ogeom_bail!(
587            NotDone,
588            "the drafted wall's fit reached {} against a target of {fit_target}",
589            fitted.error
590        );
591    }
592    Ok(fitted.curve.into())
593}
594
595/// How many stations a general draft's hinge is sampled at, all round.
596const HINGE_STATIONS: usize = 256;
597
598/// `surface` with its window run on to where it meets `neutral`, where it
599/// runs on at all: a drum's or a cone's rulings are straight lines in `v`,
600/// each meeting the plane at one closed-form `v`, and the window grows to
601/// hold them all (by a tenth of its span either way, so the crossing stands
602/// clear of its rows), short of a cone's apex. Anything else is its stored
603/// window, which a patch does not run past.
604fn reaching(
605    surface: &SurfaceGeometry,
606    neutral: Plane,
607    tol: Tolerances,
608) -> OgeomResult<SurfaceGeometry> {
609    let ((ua, ub), (va, vb)) = surface.domain();
610    let apex = match surface {
611        SurfaceGeometry::Cylinder(_) => None,
612        SurfaceGeometry::Cone(c) => {
613            let cone = c.cone();
614            Some(-cone.reference_radius() / cone.half_angle().tan())
615        }
616        _ => return Ok(surface.clone()),
617    };
618    let n = neutral.normal().vector();
619    let o = neutral.origin().to_vector();
620    let (mut lo, mut hi) = (va, vb);
621    for k in 0..HINGE_STATIONS {
622        #[allow(clippy::cast_precision_loss)]
623        let u = (ub - ua).mul_add(k as f64 / HINGE_STATIONS as f64, ua);
624        let foot = surface.point_at(u, 0.0, tol)?;
625        let ruling = surface.point_at(u, 1.0, tol)? - foot;
626        let along = n.dot(ruling);
627        if along.abs() <= tol.angular() {
628            continue;
629        }
630        let v = n.dot(o - foot.to_vector()) / along;
631        (lo, hi) = (lo.min(v), hi.max(v));
632    }
633    let margin = (hi - lo) * 0.1;
634    let (lo, hi) = (lo - margin, hi + margin);
635    // A cone runs on only to its apex: past it is the other nappe.
636    if let Some(apex) = apex
637        && ((va >= apex && lo < apex) || (vb <= apex && hi > apex))
638    {
639        return Ok(surface.clone());
640    }
641    Ok(match surface {
642        SurfaceGeometry::Cylinder(c) => {
643            ogeom_geom::CylinderSurface::new(c.cylinder(), (lo, hi))?.into()
644        }
645        SurfaceGeometry::Cone(c) => ogeom_geom::ConeSurface::new(c.cone(), (lo, hi))?.into(),
646        _ => surface.clone(),
647    })
648}
649
650/// The segments a plane cuts from a face's mesh, one per triangle it
651/// crosses, each end in the face's chart and in space. A vertex the plane
652/// passes through (the hinge running along a rim, as a draft about a base
653/// does) is a crossing in itself, not a sign to read; read as one, rounding
654/// gives every rim triangle a hair of a segment pointing anywhere. `on` is
655/// how near the plane a point counts as on it.
656fn crossing_segments(
657    mesh: &ogeom_topo::Triangulation,
658    neutral: Plane,
659    on: f64,
660) -> Vec<[((f64, f64), Point); 2]> {
661    let side: Vec<f64> = mesh
662        .positions
663        .iter()
664        .map(|p| neutral.signed_distance_to(*p))
665        .collect();
666    let mut segments: Vec<[((f64, f64), Point); 2]> = Vec::new();
667    for t in &mesh.triangles {
668        let mut ends: Vec<((f64, f64), Point)> = Vec::with_capacity(2);
669        for &corner in t {
670            let i = corner as usize;
671            if side[i].abs() <= on {
672                ends.push((mesh.parameters[i], mesh.positions[i]));
673            }
674        }
675        for k in 0..3 {
676            let (i, j) = (t[k] as usize, t[(k + 1) % 3] as usize);
677            let (a, b) = (side[i], side[j]);
678            if a.abs() <= on || b.abs() <= on || (a < 0.0) == (b < 0.0) {
679                continue;
680            }
681            let f = a / (a - b);
682            let (pa, pb) = (mesh.parameters[i], mesh.parameters[j]);
683            let (qa, qb) = (mesh.positions[i], mesh.positions[j]);
684            ends.push((
685                (pa.0 + (pb.0 - pa.0) * f, pa.1 + (pb.1 - pa.1) * f),
686                qa + (qb - qa) * f,
687            ));
688        }
689        ends.dedup_by(|a, b| a.1.distance(b.1) <= on);
690        if ends.len() == 2 && ends[0].1.distance(ends[1].1) > on {
691            segments.push([ends[0], ends[1]]);
692        }
693    }
694    segments
695}
696
697/// The turned support for any face: the ruled surface through the face's
698/// crossing with the neutral plane, its rulings the pull direction turned
699/// about the crossing's tangent by the draft: what a mould-maker means by
700/// a draft, and what the planar and revolved paths are the closed forms
701/// of. The crossing is read off the face's own mesh and corrected onto the
702/// surface, so it lies inside the face whatever the surface's chart does;
703/// a face the plane misses is read off its surface beyond it instead (see
704/// [`reaching`]), as the closed forms read their extensions.
705#[allow(clippy::too_many_arguments, reason = "one construction, all its data")]
706fn general_draft(
707    model: &Model,
708    face: &Shape,
709    surface: &SurfaceGeometry,
710    sign: f64,
711    neutral: Plane,
712    pull: Direction,
713    angle: f64,
714    tol: Tolerances,
715) -> OgeomResult<SurfaceGeometry> {
716    let n = neutral.normal().vector();
717    let mesh = ogeom_mesh::triangulate_face(
718        model,
719        face,
720        ogeom_mesh::Deflection {
721            chord: 0.05,
722            angular: 0.2,
723            ..ogeom_mesh::Deflection::default()
724        },
725        tol,
726    )?;
727    let extent = {
728        let b = mesh
729            .positions
730            .iter()
731            .fold(ogeom_math::Aabb::EMPTY, |acc, p| acc.with_point(*p));
732        match (b.low(), b.high()) {
733            (Some(lo), Some(hi)) => (hi - lo).magnitude(),
734            _ => ogeom_bail!(Construction, "the drafted face has no extent"),
735        }
736    };
737    if extent <= tol.confusion() {
738        ogeom_bail!(Construction, "the drafted face has no extent");
739    }
740
741    // The crossing, one segment per triangle the plane cuts, in the chart
742    // with its ends in space. A face the plane misses is drafted about the
743    // plane's crossing with its surface beyond it, as a planar face turns
744    // about the line its plane's extension meets the neutral one in: the
745    // same surface, meshed whole on a face of its own.
746    let on = tol.confusion() * 10.0;
747    let mut segments = crossing_segments(&mesh, neutral, on);
748    let beyond = segments.is_empty();
749    // The surface the hinge is read on: the face's own, or where the face
750    // misses the plane, the same surface run on to meet it.
751    let reached;
752    let surface: &SurfaceGeometry = if beyond {
753        reached = reaching(surface, neutral, tol)?;
754        &reached
755    } else {
756        surface
757    };
758    if beyond {
759        let mut scratch = Model::new();
760        let whole = ogeom_algo::make_natural_face(&mut scratch, surface.clone())?.shape;
761        let whole_mesh = ogeom_mesh::triangulate_face(
762            &scratch,
763            &whole,
764            ogeom_mesh::Deflection {
765                chord: 0.05,
766                angular: 0.2,
767                ..ogeom_mesh::Deflection::default()
768            },
769            tol,
770        )?;
771        segments = crossing_segments(&whole_mesh, neutral, on);
772    }
773    if segments.is_empty() {
774        ogeom_bail!(
775            Construction,
776            "the neutral plane meets neither the drafted face nor its surface; \
777             there is no hinge to turn about"
778        );
779    }
780    // Chained end to end into one run, closed or open, by where the ends
781    // are in space; across a seam the chart says two things and space one.
782    // Two runs is a face the plane crosses twice, which has no one hinge.
783    // A rim's segments come once per triangle on either side of it, so a
784    // segment already covered by the chain is dropped rather than chained.
785    let same = |a: Point, b: Point| a.distance(b) <= on;
786    let mut chain: Vec<((f64, f64), Point)> = vec![segments[0][0], segments[0][1]];
787    let mut used = vec![false; segments.len()];
788    used[0] = true;
789    loop {
790        let mut grew = false;
791        for (k, seg) in segments.iter().enumerate() {
792            if used[k] {
793                continue;
794            }
795            // The ends as they stand now: a pass grows the chain segment by
796            // segment, and the one closing it may come after the segments
797            // that brought its ends to the head and tail.
798            let tail = chain[chain.len() - 1].1;
799            let head = chain[0].1;
800            // The segment that joins the tail back to the head closes the
801            // run; it is kept, and the run is closed by it.
802            if (same(seg[0].1, tail) && same(seg[1].1, head))
803                || (same(seg[1].1, tail) && same(seg[0].1, head))
804            {
805                chain.push(chain[0]);
806                used[k] = true;
807                grew = true;
808                break;
809            }
810            let covered = |p: Point| chain.iter().any(|c| same(c.1, p));
811            if covered(seg[0].1) && covered(seg[1].1) {
812                used[k] = true;
813                continue;
814            }
815            if same(seg[0].1, tail) {
816                chain.push(seg[1]);
817            } else if same(seg[1].1, tail) {
818                chain.push(seg[0]);
819            } else if same(seg[0].1, head) {
820                chain.insert(0, seg[1]);
821            } else if same(seg[1].1, head) {
822                chain.insert(0, seg[0]);
823            } else {
824                continue;
825            }
826            used[k] = true;
827            grew = true;
828        }
829        if !grew {
830            break;
831        }
832    }
833    if used.iter().any(|u| !u) {
834        // Each run would hinge the part of the face beside it, and a part
835        // lying beside two runs would turn about both: the face would split
836        // where the two turned supports meet, which a draft (the same faces
837        // on turned supports) does not do.
838        let what = if beyond {
839            "the drafted face's surface"
840        } else {
841            "the drafted face"
842        };
843        ogeom_bail!(
844            Construction,
845            "the neutral plane crosses {what} more than once; there is no one \
846             hinge to turn about"
847        );
848    }
849    let closed = chain.len() > 3 && same(chain[0].1, chain[chain.len() - 1].1);
850    if closed {
851        chain.pop();
852    }
853    if chain.len() < 2 {
854        ogeom_bail!(
855            Construction,
856            "the neutral plane touches the drafted face at a point; there is \
857             no hinge to turn about"
858        );
859    }
860    let chain: Vec<(f64, f64)> = chain.into_iter().map(|c| c.0).collect();
861
862    // The mesh's stations are a chord apart; a cubic fitted through them
863    // sits a fraction of that chord off the true hinge. Resampled between
864    // them in the chart (across a periodic seam by the short way), and
865    // corrected onto the surface below, the stations are as dense as the
866    // fit's target wants.
867    let ((ua, ub), (va, vb)) = surface.domain();
868    // Closed counts as periodic here: a fitted tube closes on itself
869    // without repeating, and a chain crossing its join must still take
870    // the short way round.
871    let period = (
872        (surface.is_periodic_u() || surface.is_closed_u(tol)).then_some(ub - ua),
873        (surface.is_periodic_v() || surface.is_closed_v(tol)).then_some(vb - va),
874    );
875    let short = |a: f64, b: f64, period: Option<f64>| -> f64 {
876        let d = b - a;
877        match period {
878            Some(p) if d.abs() > p * 0.5 => d - p * d.signum(),
879            _ => d,
880        }
881    };
882    let chain: Vec<(f64, f64)> = {
883        // A closed hinge starts where it crosses the chart's own seam
884        // column, so the drafted support's seam stands where the old one
885        // did and the rebuild finds it there.
886        let chain: Vec<(f64, f64)> = if closed {
887            // The station exactly on the column: where the chain's segments
888            // cross `u = ua`, interpolated there; the nearest chain point
889            // otherwise.
890            let n = chain.len();
891            let mut exact: Option<(usize, (f64, f64))> = None;
892            for i in 0..n {
893                let (a, b) = (chain[i], chain[(i + 1) % n]);
894                // `b` is read as a step on from `a`: read on its own, the
895                // short way from the column flips sign half a period round
896                // as well, where the chain does not cross the column.
897                let da = short(ua, a.0, period.0);
898                let db = da + short(a.0, b.0, period.0);
899                if da == 0.0 {
900                    exact = Some((i, a));
901                    break;
902                }
903                if (da < 0.0) != (db < 0.0) && (da - db).abs() > 0.0 {
904                    let f = da / (da - db);
905                    let dv = short(a.1, b.1, period.1);
906                    exact = Some((i + 1, (ua, a.1 + dv * f)));
907                    break;
908                }
909            }
910            let (start, inserted) = exact.unwrap_or_else(|| {
911                let mut best = (0usize, f64::INFINITY);
912                for (i, c) in chain.iter().enumerate() {
913                    let d = short(ua, c.0, period.0).abs();
914                    if d < best.1 {
915                        best = (i, d);
916                    }
917                }
918                (best.0, chain[best.0])
919            });
920            let mut rotated: Vec<(f64, f64)> = Vec::with_capacity(n + 1);
921            rotated.push(inserted);
922            for k in 0..n {
923                let c = chain[(start + k) % n];
924                if rotated.len() == 1 && c == inserted {
925                    continue;
926                }
927                rotated.push(c);
928            }
929            rotated
930        } else {
931            chain
932        };
933        // Evenly by arc length, so the fit's parameter is the hinge's own
934        // length: the mesh's segments run from a hair to a chord, and a
935        // cubic through stations that uneven wanders between them.
936        let pairs = if closed { chain.len() } else { chain.len() - 1 };
937        let mut lengths = Vec::with_capacity(pairs);
938        let mut total = 0.0;
939        for i in 0..pairs {
940            let (a, b) = (chain[i], chain[(i + 1) % chain.len()]);
941            let step = surface
942                .point_at(a.0, a.1, tol)?
943                .distance(surface.point_at(b.0, b.1, tol)?);
944            lengths.push(step);
945            total += step;
946        }
947        let count = if closed {
948            HINGE_STATIONS
949        } else {
950            HINGE_STATIONS + 1
951        };
952        let mut dense = Vec::with_capacity(count);
953        let (mut pair, mut walked) = (0usize, 0.0_f64);
954        for k in 0..count {
955            #[allow(clippy::cast_precision_loss)]
956            let target = total * k as f64 / HINGE_STATIONS as f64;
957            while pair + 1 < pairs && walked + lengths[pair] < target {
958                walked += lengths[pair];
959                pair += 1;
960            }
961            let (a, b) = (chain[pair], chain[(pair + 1) % chain.len()]);
962            let (du, dv) = (short(a.0, b.0, period.0), short(a.1, b.1, period.1));
963            let f = if lengths[pair] > 0.0 {
964                ((target - walked) / lengths[pair]).clamp(0.0, 1.0)
965            } else {
966                0.0
967            };
968            dense.push((a.0 + du * f, a.1 + dv * f));
969        }
970        dense
971    };
972
973    // Each station corrected onto the surface's crossing with the plane
974    // (the mesh's chord is not the surface), and read for its tangent and
975    // outward normal.
976    // A chart position corrected onto the crossing (along `v` alone where
977    // `pinned`), with the point there, the crossing's direction (either
978    // way along it) and the outward normal.
979    let on_hinge = |(mut u, mut v): (f64, f64),
980                    pinned: bool|
981     -> OgeomResult<((f64, f64), Point, Vector, Vector)> {
982        for _ in 0..8 {
983            let p = surface.point_at(u, v, tol)?;
984            let f = neutral.signed_distance_to(p);
985            if f.abs() <= tol.confusion() * 1e-2 {
986                break;
987            }
988            let (du, dv) = surface.d1_at(u, v, tol)?;
989            let g = (if pinned { 0.0 } else { n.dot(du) }, n.dot(dv));
990            let g2 = g.0 * g.0 + g.1 * g.1;
991            if g2 <= 0.0 {
992                break;
993            }
994            u -= f * g.0 / g2;
995            v -= f * g.1 / g2;
996        }
997        let p = surface.point_at(u, v, tol)?;
998        let (du, dv) = surface.d1_at(u, v, tol)?;
999        let raw = du.cross(dv);
1000        if raw.magnitude() <= tol.confusion() {
1001            ogeom_bail!(Construction, "the drafted face has no normal on its hinge");
1002        }
1003        let outward = raw / raw.magnitude() * sign;
1004        // Along the crossing: square to both normals.
1005        let t = outward.cross(n);
1006        if t.magnitude() <= tol.angular() {
1007            ogeom_bail!(
1008                Construction,
1009                "the neutral plane is tangent to the drafted face; there is no \
1010                 hinge to turn about"
1011            );
1012        }
1013        Ok(((u, v), p, t / t.magnitude(), outward))
1014    };
1015    let mut feet: Vec<(f64, f64)> = Vec::with_capacity(chain.len());
1016    let mut hinges: Vec<Point> = Vec::with_capacity(chain.len());
1017    let mut tangents: Vec<Vector> = Vec::with_capacity(chain.len());
1018    let mut outwards: Vec<Vector> = Vec::with_capacity(chain.len());
1019    for (k, &at) in chain.iter().enumerate() {
1020        // The first station of a closed hinge is the seam column's, and
1021        // stays on it.
1022        let (foot, p, mut t, outward) = on_hinge(at, closed && k == 0)?;
1023        // Run the way the chain runs.
1024        let next = chain[(k + 1) % chain.len()];
1025        let prev = chain[(k + chain.len() - 1) % chain.len()];
1026        let ahead =
1027            surface.point_at(next.0, next.1, tol)? - surface.point_at(prev.0, prev.1, tol)?;
1028        if t.dot(ahead) < 0.0 {
1029            t = -t;
1030        }
1031        feet.push(foot);
1032        hinges.push(p);
1033        tangents.push(t);
1034        outwards.push(outward);
1035    }
1036
1037    // The sense, probed at the middle station as the other paths probe:
1038    // the turn whose outward normal leans furthest towards the pull is the
1039    // inward one, and the caller's sign picks inward or outward through it.
1040    let m = hinges.len() / 2;
1041    let axis_m = ogeom_math::Axis::new(hinges[m], Direction::new(tangents[m], tol)?);
1042    let mut leaning = 1.0;
1043    let mut best = f64::NEG_INFINITY;
1044    for sense in [1.0_f64, -1.0] {
1045        let turn = Transform::rotation(axis_m, angle.abs() * sense);
1046        let lean = turn.apply_vector(outwards[m]).dot(pull.vector());
1047        if lean > best {
1048            best = lean;
1049            leaning = sense;
1050        }
1051    }
1052    let theta = angle * leaning;
1053
1054    // One ruling per station: the pull turned about the hinge's tangent.
1055    let mut rulings: Vec<Vector> = Vec::with_capacity(hinges.len());
1056    for (hinge, tangent) in hinges.iter().zip(&tangents) {
1057        let turn = Transform::rotation(
1058            ogeom_math::Axis::new(*hinge, Direction::new(*tangent, tol)?),
1059            theta,
1060        );
1061        rulings.push(turn.apply_vector(pull.vector()));
1062    }
1063    // How far the face reaches along the pull either side of its hinge,
1064    // grown with the turn so the wall still reaches the neighbours it
1065    // re-meets.
1066    // Measured from every hinge station, not the middle one: an oblique
1067    // hinge rises and falls along the pull, and a ruling from its low end
1068    // must still reach the face's top.
1069    let (mut s_lo, mut s_hi) = (f64::INFINITY, f64::NEG_INFINITY);
1070    let (mut h_lo, mut h_hi) = (f64::INFINITY, f64::NEG_INFINITY);
1071    for h in &hinges {
1072        let s = h.to_vector().dot(pull.vector());
1073        h_lo = h_lo.min(s);
1074        h_hi = h_hi.max(s);
1075    }
1076    for p in &mesh.positions {
1077        let s = p.to_vector().dot(pull.vector());
1078        s_lo = s_lo.min(s - h_hi);
1079        s_hi = s_hi.max(s - h_lo);
1080    }
1081    let grow = extent.mul_add(0.5, 1.0) * angle.abs().tan() + tol.confusion();
1082    let (s_lo, s_hi) = (s_lo - grow, s_hi + grow);
1083    // How many stations continue an open hinge before its first.
1084    let mut lead = 0;
1085    if !closed {
1086        // An open hinge is continued straight past both ends, for the
1087        // same reason.
1088        let extend = |hinges: &mut Vec<Point>, rulings: &mut Vec<Vector>, front: bool| -> usize {
1089            let (i0, i1) = if front {
1090                (0, 1)
1091            } else {
1092                (hinges.len() - 1, hinges.len() - 2)
1093            };
1094            let d1 = hinges[i0] - hinges[i1];
1095            let steps = (grow / d1.magnitude().max(tol.confusion())).ceil().max(2.0);
1096            #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
1097            let steps = (steps as usize).min(16);
1098            for k in 1..=steps {
1099                #[allow(clippy::cast_precision_loss)]
1100                let station = (hinges[i0] + d1 * k as f64, rulings[i0]);
1101                if front {
1102                    hinges.insert(0, station.0);
1103                    rulings.insert(0, station.1);
1104                } else {
1105                    hinges.push(station.0);
1106                    rulings.push(station.1);
1107                }
1108            }
1109            steps
1110        };
1111        lead = extend(&mut hinges, &mut rulings, true);
1112        extend(&mut hinges, &mut rulings, false);
1113    } else {
1114        hinges.push(hinges[0]);
1115        rulings.push(rulings[0]);
1116    }
1117    for edge in [s_lo, s_hi] {
1118        for i in 0..hinges.len() - 1 {
1119            let step = (hinges[i + 1] + rulings[i + 1] * edge) - (hinges[i] + rulings[i] * edge);
1120            if step.dot(hinges[i + 1] - hinges[i]) <= 0.0 {
1121                ogeom_bail!(
1122                    Construction,
1123                    "the draft folds the wall onto itself inside the drafted \
1124                     window; refused; see docs/CAPABILITIES.md, offset.draft"
1125                );
1126            }
1127        }
1128    }
1129    // The ruled surface itself, exactly: the wall's two border rows, at the
1130    // window's two heights along the rulings, fitted as cubics on one knot
1131    // vector at the *same* parameters and checked between the stations,
1132    // and the surface linear between them: degree one along the ruling, so
1133    // a straight line is a straight line, and every point of the window
1134    // stands between two rows each within the target. A grid fit through
1135    // rows at several heights parameterizes each row by its own chord and
1136    // averages, and rows that converge along their rulings disagree by
1137    // enough for a cubic across them to wander.
1138    // Chord-length parameters, not centripetal: the stations are as evenly
1139    // spaced as the mesh's segments let them be, and centripetal
1140    // parameters kink wherever the spacing changes, which a cubic then
1141    // cannot follow. By chord the parameter is the arc length whatever
1142    // the spacing.
1143    let params: Vec<f64> = {
1144        let mut out = Vec::with_capacity(hinges.len());
1145        let mut total = 0.0;
1146        out.push(0.0);
1147        for pair in hinges.windows(2) {
1148            total += pair[0].distance(pair[1]);
1149            out.push(total);
1150        }
1151        if total > 0.0 {
1152            for t in &mut out {
1153                *t /= total;
1154            }
1155        }
1156        if let Some(last) = out.last_mut() {
1157            *last = 1.0;
1158        }
1159        out
1160    };
1161    // The hinge and its ruling anywhere along it: a station's own, and
1162    // between two stations the crossing found again from the chart between
1163    // their feet, its ruling turned there as at a station; along an open
1164    // hinge's straight continuation, the line and its end's ruling.
1165    let stations = feet.len();
1166    let wrap = |x: f64, lo: f64, period: Option<f64>| -> f64 {
1167        period.map_or(x, |p| lo + (x - lo).rem_euclid(p))
1168    };
1169    let exact_at = |at: f64| -> OgeomResult<(Point, Vector)> {
1170        let i = params
1171            .partition_point(|p| *p <= at)
1172            .saturating_sub(1)
1173            .min(params.len() - 2);
1174        let f = (at - params[i]) / (params[i + 1] - params[i]);
1175        if f <= 0.0 {
1176            return Ok((hinges[i], rulings[i]));
1177        }
1178        let (k0, k1) = (i.wrapping_sub(lead), (i + 1).wrapping_sub(lead));
1179        if k0 >= stations || k1 > stations || (!closed && k1 == stations) {
1180            return Ok((hinges[i] + (hinges[i + 1] - hinges[i]) * f, rulings[i]));
1181        }
1182        let (a, b) = (feet[k0], feet[k1 % stations]);
1183        let (mut u, mut v) = (
1184            wrap(a.0 + short(a.0, b.0, period.0) * f, ua, period.0),
1185            wrap(a.1 + short(a.1, b.1, period.1) * f, va, period.1),
1186        );
1187        // On the plane, and as far along the chord between the two
1188        // stations as the parameter is between theirs: a correction onto
1189        // the plane alone would slide along the hinge as it went, and the
1190        // hinge's pace would kink at every station.
1191        let (from, chord) = (hinges[i], hinges[i + 1] - hinges[i]);
1192        for _ in 0..8 {
1193            let q = surface.point_at(u, v, tol)?;
1194            let g = (
1195                neutral.signed_distance_to(q),
1196                (q - from).dot(chord) - f * chord.dot(chord),
1197            );
1198            if g.0.abs() <= tol.confusion() * 1e-2
1199                && g.1.abs() <= tol.confusion() * 1e-2 * chord.magnitude()
1200            {
1201                break;
1202            }
1203            let (du, dv) = surface.d1_at(u, v, tol)?;
1204            let (a11, a12, a21, a22) = (n.dot(du), n.dot(dv), chord.dot(du), chord.dot(dv));
1205            let det = a11 * a22 - a12 * a21;
1206            if det.abs() <= f64::MIN_POSITIVE {
1207                break;
1208            }
1209            u = wrap(u - (g.0 * a22 - a12 * g.1) / det, ua, period.0);
1210            v = wrap(v - (a11 * g.1 - a21 * g.0) / det, va, period.1);
1211        }
1212        let (_, p, t, _) = on_hinge((u, v), false)?;
1213        let t = if t.dot(tangents[k0]) < 0.0 { -t } else { t };
1214        let turn = Transform::rotation(ogeom_math::Axis::new(p, Direction::new(t, tol)?), theta);
1215        Ok((p, turn.apply_vector(pull.vector())))
1216    };
1217    // A closed hinge leaves its first station and comes back to it on the
1218    // two sides of the face's seam, where a face closed only to its
1219    // position turns its normal, and its rulings with it. The wall closes
1220    // on one ruling there, the mean of the two sides, and turns smoothly
1221    // from it to each side's own over a thirty-second of the way round
1222    // either side: a turn within one station would ask the fit for a knot
1223    // at every sample there. Each side continued past the seam on its
1224    // tangent plane meets the other near the mean ruling (on it, to first
1225    // order, where the hinge's tangent is square to the pull), so the
1226    // blend keeps the wall on the drafted surface; it is each side's last
1227    // ruling the wall leaves, by up to half the turn times the reach.
1228    const SEAM_BLEND: f64 = 1.0 / 32.0;
1229    let seam = if closed {
1230        let leave = exact_at(params[1] * 1e-9)?.1;
1231        let back = exact_at(1.0 - (1.0 - params[params.len() - 2]) * 1e-9)?.1;
1232        let mean = leave + back;
1233        Some((leave, back, mean / mean.magnitude()))
1234    } else {
1235        None
1236    };
1237    let ruled_at = |at: f64| -> OgeomResult<(Point, Vector)> {
1238        let Some((leave, back, mean)) = seam else {
1239            return exact_at(at);
1240        };
1241        if at <= 0.0 || at >= 1.0 {
1242            return Ok((hinges[0], mean));
1243        }
1244        let (p, r) = exact_at(at)?;
1245        // The share of the turn left: one at the seam, none past the blend,
1246        // flat at both ends.
1247        let left = |x: f64| {
1248            let x = x.clamp(0.0, 1.0);
1249            1.0 - x * x * 2.0f64.mul_add(-x, 3.0)
1250        };
1251        let r = if at < SEAM_BLEND {
1252            r + (mean - leave) * left(at / SEAM_BLEND)
1253        } else if at > 1.0 - SEAM_BLEND {
1254            r + (mean - back) * left((1.0 - at) / SEAM_BLEND)
1255        } else {
1256            return Ok((p, r));
1257        };
1258        Ok((p, r / r.magnitude()))
1259    };
1260    // Both rows on one knot vector, refined where either strays between
1261    // the samples: a surface fit through the two rows, which is linear
1262    // across them.
1263    let fit_target = (tol.confusion() * 1e3).max(1e-4);
1264    let rows = ogeom_geom::fit::fit_surface_sampled(
1265        |at, side| {
1266            let (h, r) = ruled_at(at)?;
1267            Ok(h + r * if side < 0.5 { s_lo } else { s_hi })
1268        },
1269        &ogeom_geom::fit::Sampling {
1270            us: params.clone(),
1271            vs: vec![0.0, 1.0],
1272            between: (true, false),
1273            closed_v: false,
1274            most: 1024,
1275        },
1276        3,
1277        fit_target,
1278        tol,
1279    )?;
1280    if !rows.met {
1281        ogeom_bail!(
1282            NotDone,
1283            "the drafted wall's border rows reached {} against a target of {fit_target}",
1284            rows.error
1285        );
1286    }
1287    let fitted = rows.curve;
1288    let (k, l) = (fitted.grid().u_count(), fitted.grid().v_count());
1289    if l != 2 || fitted.v_knots().degree() != 1 {
1290        ogeom_bail!(
1291            Construction,
1292            "the drafted wall's border rows did not fit as a ruled surface"
1293        );
1294    }
1295    let corners: Vec<Point> = fitted.grid().points().iter().map(|w| w.point()).collect();
1296    let (lc, hc): (Vec<Point>, Vec<Point>) =
1297        (0..k).map(|i| (corners[i * 2], corners[i * 2 + 1])).unzip();
1298    // The face keeps its orientation flag, so the new surface's own normal
1299    // must turn the way the old one's did; whether it does depends on which
1300    // way the hinge chained. Measured at the hinge's middle, where the two
1301    // surfaces meet, and the rulings run from the far end back if not.
1302    let mid = {
1303        let (ud, _) = fitted.domain();
1304        f64::midpoint(ud.0, ud.1)
1305    };
1306    let (low_mid, high_mid) = (
1307        fitted.point_at(mid, 0.0, tol)?,
1308        fitted.point_at(mid, 1.0, tol)?,
1309    );
1310    let hinge_mid = low_mid + (high_mid - low_mid) * (-s_lo / (s_hi - s_lo));
1311    let across = fitted.d1_at(mid, 0.0, tol)?.0;
1312    let old_normal = {
1313        let foot = ogeom_algo::project_on_surface(surface, hinge_mid, 16, tol)?;
1314        let (u, v) = foot.parameters;
1315        surface.normal_at(u, v, tol)?.vector()
1316    };
1317    // du x dv of the new surface at the hinge: along the hinge, crossed
1318    // with up the ruling.
1319    let agrees = across.cross(high_mid - low_mid).dot(old_normal) >= 0.0;
1320    let (first, second, v_range) = if agrees {
1321        (&lc, &hc, (s_lo, s_hi))
1322    } else {
1323        (&hc, &lc, (-s_hi, -s_lo))
1324    };
1325    let mut net: Vec<Point> = Vec::with_capacity(k * 2);
1326    for (a, b) in first.iter().zip(second) {
1327        net.push(*a);
1328        net.push(*b);
1329    }
1330    let grid = ogeom_math::ControlGrid::new(net, k, 2)?;
1331    // The chart: `u` over the old surface's own `u` domain for a closed
1332    // hinge, so the seam column carries over; `v` the height along the
1333    // ruling in the model's own units.
1334    let u_knots = if closed {
1335        fitted.u_knots().reparameterized(ua, ub)?
1336    } else {
1337        fitted.u_knots().clone()
1338    };
1339    let v_knots =
1340        ogeom_math::KnotVector::clamped_uniform(1, 2)?.reparameterized(v_range.0, v_range.1)?;
1341    Ok(ogeom_geom::BSplineSurface::new(u_knots, v_knots, &grid, tol)?.into())
1342}
1343
1344/// Which sign turns a surface's raw normal (du x dv) outward, read from
1345/// the solid itself: probed a step off the face midpoint on both sides, at
1346/// growing steps until one side is material and the other is not.
1347///
1348/// # Errors
1349///
1350/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if no
1351/// probe separates the sides: a wall thinner than the probe can resolve.
1352fn outward_sign(
1353    model: &Model,
1354    solid: &Shape,
1355    face: &Shape,
1356    surface: &SurfaceGeometry,
1357    tol: Tolerances,
1358) -> OgeomResult<f64> {
1359    use ogeom_algo::Containment;
1360    // A point genuinely on the face (the surface's domain midpoint may lie
1361    // outside the trim), from the face's own triangulation, at its largest
1362    // triangle's centre.
1363    let mesh = ogeom_mesh::triangulate_face(model, face, ogeom_mesh::Deflection::default(), tol)?;
1364    let mut at = None;
1365    let mut largest = 0.0_f64;
1366    for t in &mesh.triangles {
1367        let [a, b, c] = [
1368            mesh.positions[t[0] as usize],
1369            mesh.positions[t[1] as usize],
1370            mesh.positions[t[2] as usize],
1371        ];
1372        let area = (b - a).cross(c - a).magnitude();
1373        if area > largest {
1374            largest = area;
1375            let params = [
1376                mesh.parameters[t[0] as usize],
1377                mesh.parameters[t[1] as usize],
1378                mesh.parameters[t[2] as usize],
1379            ];
1380            at = Some((
1381                (params[0].0 + params[1].0 + params[2].0) / 3.0,
1382                (params[0].1 + params[1].1 + params[2].1) / 3.0,
1383            ));
1384        }
1385    }
1386    let Some((um, vm)) = at else {
1387        ogeom_bail!(Construction, "the drafted face has no interior to probe");
1388    };
1389    let p = surface.point_at(um, vm, tol)?;
1390    let (du, dv) = surface.d1_at(um, vm, tol)?;
1391    let n = du.cross(dv);
1392    let m = n.magnitude();
1393    if m <= tol.confusion() {
1394        ogeom_bail!(Construction, "the face has no normal at its midpoint");
1395    }
1396    let n = n / m;
1397    let scale = largest.sqrt().max(tol.confusion() * 1e3);
1398    for eps_scale in [1e-3, 1e-2, 5e-2] {
1399        let eps = scale * eps_scale;
1400        let deflection = ogeom_mesh::Deflection {
1401            chord: (eps * 0.1).max(1e-4),
1402            ..ogeom_mesh::Deflection::default()
1403        };
1404        let ahead = ogeom_algo::classify_in_solid(model, solid, p + n * eps, deflection, tol)?;
1405        let behind = ogeom_algo::classify_in_solid(model, solid, p - n * eps, deflection, tol)?;
1406        match (ahead, behind) {
1407            (Containment::Out, Containment::In) => return Ok(1.0),
1408            (Containment::In, Containment::Out) => return Ok(-1.0),
1409            _ => {}
1410        }
1411    }
1412    ogeom_bail!(
1413        Construction,
1414        "cannot read which side of the drafted face holds material; the          wall is thinner than the probe can resolve"
1415    )
1416}
1417
1418/// A point on the line where two planes meet, nearest their origins.
1419fn meet(a: Plane, b: Plane, along: Vector, tol: Tolerances) -> OgeomResult<Point> {
1420    let rows = [a.normal().vector(), b.normal().vector(), along];
1421    let rhs = [
1422        rows[0].dot(a.origin().to_vector()),
1423        rows[1].dot(b.origin().to_vector()),
1424        along.dot(Point::midpoint(a.origin(), b.origin()).to_vector()),
1425    ];
1426    let det = rows[0].dot(rows[1].cross(rows[2]));
1427    if det.abs() <= tol.confusion() {
1428        ogeom_bail!(Construction, "the two planes do not meet in a line");
1429    }
1430    Ok(Point::ORIGIN
1431        + (rows[1].cross(rows[2]) * rhs[0]
1432            + rows[2].cross(rows[0]) * rhs[1]
1433            + rows[0].cross(rows[1]) * rhs[2])
1434            / det)
1435}