Skip to main content

ogeom_offset/
sweep.rs

1//! Pipes and lofts: the sweeps whose surfaces are already in the vocabulary.
2//!
3//! A circular profile along a straight spine is a cylinder; around a full
4//! circle it is a torus; along an arc it is a torus segment, built from two
5//! half-tube patches and two meridian caps so that every edge is a circle
6//! with a closed-form chart. A ruled loft between two parallel sections is
7//! walls of planes and cones: segment to segment gives the planar quad,
8//! coaxial circle to circle gives the frustum the cone primitive already
9//! builds. The sweeps that need *new* surfaces (free-form spines, skew
10//! ruled walls, smoothed skinning through many sections) are recorded in
11//! docs/PARITY.md (offset.sweeps), not approximated here.
12
13use ogeom_algo::{
14    Built, History, edge_vertices, make_cone, make_cylinder, make_edge, make_edge_between,
15    make_face_with_pcurves, make_solid, make_torus, make_vertex, sew,
16};
17use ogeom_core::{OgeomResult, Tolerances, ogeom_bail};
18use ogeom_geom::Curve3d as _;
19use ogeom_geom::Transformable as _;
20use ogeom_geom::{
21    CircleCurve, Curve, Line2d, LineCurve, PlaneSurface, SurfaceGeometry, TorusSurface,
22};
23use ogeom_math::{Circle, Direction, Frame, Plane, Point, Point2, Torus, Transform, Vector};
24use ogeom_topo::{EdgeData, EdgeRepr, Filter, Model, Shape, ShapeType, VertexData, explore};
25
26/// Sweep a circular profile of `radius` along a spine edge.
27///
28/// A straight spine gives a cylinder, a full circular spine a torus, an arc
29/// a torus segment. The history generates the solid from the spine.
30///
31/// # Errors
32///
33/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the spine is
34/// not straight or circular, the profile radius is not usable, or the tube
35/// would swallow its own spine.
36pub fn make_pipe(
37    model: &mut Model,
38    spine: &Shape,
39    radius: f64,
40    tol: Tolerances,
41) -> OgeomResult<Built> {
42    if !radius.is_finite() || radius <= tol.confusion() {
43        ogeom_bail!(Construction, "a pipe of radius {radius} holds nothing");
44    }
45    let (curve, range) = {
46        let Some(data) = model.node(spine).and_then(|n| n.data().as_edge()) else {
47            ogeom_bail!(Construction, "a pipe runs along an edge");
48        };
49        let Some(EdgeRepr::Curve3d { curve, range, .. }) = data.curve3d() else {
50            ogeom_bail!(Construction, "the spine has no curve");
51        };
52        let Some(geometry) = model.geometry().curve(*curve) else {
53            ogeom_bail!(Dangling, "curve is not in this model");
54        };
55        (geometry.clone(), *range)
56    };
57    let mut built = match &curve {
58        Curve::Line(line) => {
59            let start = curve.point_at(range.0, tol)?;
60            let length = range.1 - range.0;
61            let frame = Frame::about(start, line.axis().direction);
62            make_cylinder(model, frame, radius, length, tol)?
63        }
64        Curve::Circle(c) => {
65            let circle = c.circle();
66            if radius >= circle.radius() - tol.confusion() {
67                ogeom_bail!(
68                    Construction,
69                    "a tube of radius {radius} swallows its spine of radius {}",
70                    circle.radius()
71                );
72            }
73            let closed = curve
74                .point_at(range.0, tol)?
75                .distance(curve.point_at(range.1, tol)?)
76                <= tol.confusion();
77            if closed {
78                make_torus(model, circle.frame(), circle.radius(), radius, tol)?
79            } else {
80                pipe_segment(model, circle, range, radius, tol)?
81            }
82        }
83        _ => ogeom_bail!(
84            Construction,
85            "a pipe along a free-form spine needs the sweep-surface \
86             machinery; see docs/PARITY.md, offset.sweeps"
87        ),
88    };
89    built.history.generate(spine, built.shape.clone());
90    Ok(built)
91}
92
93/// The torus segment: two half-tube patches and two meridian caps.
94///
95/// The tube circles are framed so their own parameter *is* the torus tube
96/// angle, which makes every tube pcurve a vertical line in the chart and the
97/// two patches the clean rectangles `v ∈ [0, π]` and `[π, 2π]`. The outer
98/// equator is then the seam between the halves across the period (one edge,
99/// two chart rows), which is exactly what [`ogeom_algo::attach_seam`] exists to
100/// say.
101fn pipe_segment(
102    model: &mut Model,
103    spine: Circle,
104    range: (f64, f64),
105    radius: f64,
106    tol: Tolerances,
107) -> OgeomResult<Built> {
108    let frame = spine.frame();
109    let (x, y, z) = (frame.x().vector(), frame.y().vector(), frame.z().vector());
110    let major = spine.radius();
111    let radial = |u: f64| x * u.cos() + y * u.sin();
112    let tangent = |u: f64| x * -u.sin() + y * u.cos();
113    let tube_point = |u: f64, v: f64| {
114        frame.origin() + radial(u) * radius.mul_add(v.cos(), major) + z * (radius * v.sin())
115    };
116    let pi = core::f64::consts::PI;
117    let tau = core::f64::consts::TAU;
118
119    let torus: SurfaceGeometry = TorusSurface::new(Torus::new(frame, major, radius, tol)?).into();
120    let surface_id = model.geometry_mut().add_surface(torus);
121
122    // Vertices at the tube's v = 0 and v = π points of each end.
123    let ends = [range.0, range.1];
124    let mut verts: Vec<Vec<Shape>> = Vec::new();
125    for &u in &ends {
126        verts.push(vec![
127            make_vertex(model, tube_point(u, 0.0)).shape,
128            make_vertex(model, tube_point(u, pi)).shape,
129        ]);
130    }
131
132    // The tube circles at each end, split at v = 0 and v = π, framed so that
133    // the circle's parameter equals the torus tube angle: `z` against the
134    // spine tangent makes the frame's `y` the torus's own axis.
135    let mut tube_arcs: Vec<Vec<Shape>> = Vec::new();
136    for (k, &u) in ends.iter().enumerate() {
137        let centre = frame.origin() + radial(u) * major;
138        let circle = Circle::new(
139            Frame::new(
140                centre,
141                Direction::new(-tangent(u), tol)?,
142                Direction::new(radial(u), tol)?,
143                tol,
144            )?,
145            radius,
146            tol,
147        )?;
148        let curve = Curve::Circle(CircleCurve::new(circle));
149        let arcs = vec![
150            make_edge_between(
151                model,
152                curve.clone(),
153                (0.0, pi),
154                &verts[k][0],
155                &verts[k][1],
156                tol,
157            )?
158            .shape,
159            make_edge_between(model, curve, (pi, tau), &verts[k][1], &verts[k][0], tol)?.shape,
160        ];
161        // In the chart both arcs run up the column at this end's angle.
162        let column = Line2d::over(
163            ogeom_math::Axis2::new(Point2::new(u, 0.0), ogeom_math::Direction2::Y),
164            0.0,
165            tau,
166        )?;
167        ogeom_algo::attach_pcurve(
168            model,
169            &arcs[0],
170            column.into(),
171            surface_id,
172            ogeom_topo::Location::identity(),
173            (0.0, pi),
174        )?;
175        ogeom_algo::attach_pcurve(
176            model,
177            &arcs[1],
178            column.into(),
179            surface_id,
180            ogeom_topo::Location::identity(),
181            (pi, tau),
182        )?;
183        tube_arcs.push(arcs);
184    }
185
186    // The long edges: the parallels at v = 0 and v = π, parameterized by the
187    // spine's own angle.
188    let parallel = |model: &mut Model, v: f64, from: &Shape, to: &Shape| -> OgeomResult<Shape> {
189        let height = radius * v.sin();
190        let ring = radius.mul_add(v.cos(), major);
191        let circle = Circle::new(
192            Frame::new(frame.origin() + z * height, frame.z(), frame.x(), tol)?,
193            ring,
194            tol,
195        )?;
196        let curve = Curve::Circle(CircleCurve::new(circle));
197        Ok(make_edge_between(model, curve, range, from, to, tol)?.shape)
198    };
199    let row = |v: f64| -> OgeomResult<Line2d> {
200        Line2d::over(
201            ogeom_math::Axis2::new(Point2::new(0.0, v), ogeom_math::Direction2::X),
202            range.0 - 1.0,
203            range.1 + 1.0,
204        )
205    };
206    let inner = parallel(model, pi, &verts[0][1], &verts[1][1])?;
207    ogeom_algo::attach_pcurve(
208        model,
209        &inner,
210        row(pi)?.into(),
211        surface_id,
212        ogeom_topo::Location::identity(),
213        range,
214    )?;
215    // The outer equator bounds both halves across the period: v = 2π for its
216    // forward use under the upper patch, v = 0 for its reversed use under the
217    // lower: a seam, said as one.
218    let outer = parallel(model, 0.0, &verts[0][0], &verts[1][0])?;
219    ogeom_algo::attach_seam(
220        model,
221        &outer,
222        row(tau)?.into(),
223        row(0.0)?.into(),
224        surface_id,
225        ogeom_topo::Location::identity(),
226        range,
227    )?;
228
229    // The two half-tube patches, on the one registered surface.
230    let lower = {
231        let wire = ogeom_algo::make_wire(
232            model,
233            &[
234                tube_arcs[0][0].clone(),
235                inner.clone(),
236                tube_arcs[1][0].reversed(),
237                outer.reversed(),
238            ],
239            tol,
240        )?
241        .shape;
242        ogeom_algo::make_face_on(model, surface_id, std::slice::from_ref(&wire), tol)?.shape
243    };
244    let upper = {
245        let wire = ogeom_algo::make_wire(
246            model,
247            &[
248                tube_arcs[0][1].clone(),
249                outer.clone(),
250                tube_arcs[1][1].reversed(),
251                inner.reversed(),
252            ],
253            tol,
254        )?
255        .shape;
256        ogeom_algo::make_face_on(model, surface_id, std::slice::from_ref(&wire), tol)?.shape
257    };
258
259    // The meridian caps, their outward normals along the spine and away from
260    // the material between the ends.
261    let mut caps: Vec<Shape> = Vec::new();
262    for (k, &u) in ends.iter().enumerate() {
263        let outward = if k == 0 { -tangent(u) } else { tangent(u) };
264        let centre = frame.origin() + radial(u) * major;
265        let plane = Plane::through(centre, Direction::new(outward, tol)?);
266        let reach = (major + radius) * 2.0;
267        let surface: SurfaceGeometry =
268            PlaneSurface::over(plane, (-reach, reach), (-reach, reach))?.into();
269        caps.push(
270            make_face_with_pcurves(
271                model,
272                surface,
273                &[vec![tube_arcs[k][0].clone(), tube_arcs[k][1].clone()]],
274                tol,
275            )?
276            .shape,
277        );
278    }
279
280    let faces = [lower, upper, caps[0].clone(), caps[1].clone()];
281    let sewn = sew(model, &faces, tol)?;
282    if sewn.shells.len() != 1 || !ogeom_algo::is_shell_closed(model, &sewn.shells[0])? {
283        ogeom_bail!(Construction, "the pipe segment did not close");
284    }
285    make_solid(model, std::slice::from_ref(&sewn.shells[0]))
286}
287
288/// Loft two parallel closed sections into a solid, ruled.
289///
290/// Two coaxial circles give the cylinder or the cone frustum; two polygons
291/// with the same corner count give planar walls. The sections pair edge by
292/// edge in traversal order.
293///
294/// A wall between two segments that are not coplanar is the bilinear patch
295/// through its four corners (the ruled surface between them, exact), so
296/// sections may be turned against each other or differ in shape.
297///
298/// # Errors
299///
300/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the sections
301/// are not both circles or both polygons of the same count, or are not on
302/// parallel planes.
303pub fn make_loft(
304    model: &mut Model,
305    bottom: &Shape,
306    top: &Shape,
307    tol: Tolerances,
308) -> OgeomResult<Built> {
309    // A vertex for either section is the loft to a point: a cone over a
310    // circle, an exact pyramid over a polygon.
311    match (model.kind_of(bottom)?, model.kind_of(top)?) {
312        (ShapeType::Wire, ShapeType::Vertex) => {
313            return loft_to_point(model, bottom, top, tol);
314        }
315        (ShapeType::Vertex, ShapeType::Wire) => {
316            let mut built = loft_to_point(model, top, bottom, tol)?;
317            // The apex was named first: the same solid, the history the same.
318            built.history.generate(bottom, built.shape.clone());
319            return Ok(built);
320        }
321        _ => {}
322    }
323    for wire in [bottom, top] {
324        if model.kind_of(wire)? != ShapeType::Wire {
325            ogeom_bail!(Construction, "a loft runs between wires");
326        }
327        if !ogeom_algo::is_wire_closed(model, wire, tol)? {
328            ogeom_bail!(Construction, "a loft section must be closed");
329        }
330    }
331    let circle_of = |model: &Model, wire: &Shape| -> OgeomResult<Option<Circle>> {
332        let edges = explore(model, wire, Filter::OfType(ShapeType::Edge))?;
333        if edges.len() != 1 {
334            return Ok(None);
335        }
336        let Some(data) = model.node(&edges[0]).and_then(|n| n.data().as_edge()) else {
337            return Ok(None);
338        };
339        let Some(EdgeRepr::Curve3d { curve, .. }) = data.curve3d() else {
340            return Ok(None);
341        };
342        match model.geometry().curve(*curve) {
343            Some(Curve::Circle(c)) => Ok(Some(c.circle())),
344            _ => Ok(None),
345        }
346    };
347
348    if let (Some(lower), Some(upper)) = (circle_of(model, bottom)?, circle_of(model, top)?) {
349        // Coaxial circles: the revolved primitives already build these.
350        let axis = lower.frame().z().vector();
351        let rise = upper.centre() - lower.centre();
352        let height = rise.dot(axis);
353        if rise.cross(axis).magnitude() > tol.confusion() * 10.0 || height.abs() <= tol.confusion()
354        {
355            ogeom_bail!(
356                Construction,
357                "lofted circles must be coaxial on parallel planes; the \
358                 oblique loft needs the sweep machinery; see the deferred \
359                 table"
360            );
361        }
362        let frame = if height > 0.0 {
363            Frame::new(lower.centre(), lower.frame().z(), lower.frame().x(), tol)?
364        } else {
365            Frame::new(
366                lower.centre(),
367                lower.frame().z().reversed(),
368                lower.frame().x(),
369                tol,
370            )?
371        };
372        let mut built = if (lower.radius() - upper.radius()).abs() <= tol.confusion() {
373            make_cylinder(model, frame, lower.radius(), height.abs(), tol)?
374        } else {
375            make_cone(
376                model,
377                frame,
378                lower.radius(),
379                upper.radius(),
380                height.abs(),
381                tol,
382            )?
383        };
384        built.history.generate(bottom, built.shape.clone());
385        built.history.generate(top, built.shape.clone());
386        return Ok(built);
387    }
388
389    // Polygonal sections: matched corners, planar walls.
390    let corners_of = |model: &Model, wire: &Shape| -> OgeomResult<Vec<Point>> {
391        let mut out = Vec::new();
392        for edge in explore(model, wire, Filter::OfType(ShapeType::Edge))? {
393            let Some(data) = model.node(&edge).and_then(|n| n.data().as_edge()) else {
394                ogeom_bail!(Construction, "a section edge holds no data");
395            };
396            let Some(EdgeRepr::Curve3d { curve, .. }) = data.curve3d() else {
397                ogeom_bail!(Construction, "a section edge has no curve");
398            };
399            let Some(Curve::Line(_)) = model.geometry().curve(*curve) else {
400                ogeom_bail!(
401                    Construction,
402                    "a mixed or curved section needs the skinning machinery; see \
403                     docs/PARITY.md, offset.sweeps"
404                );
405            };
406            let Some((sv, _)) = edge_vertices(model, &edge)? else {
407                ogeom_bail!(Construction, "a section edge has no vertices");
408            };
409            let Some(data) = model.node(&sv).and_then(|n| n.data().as_vertex()) else {
410                ogeom_bail!(Construction, "a section vertex holds no point");
411            };
412            out.push(sv.transform(model.datums())?.apply(data.point));
413        }
414        Ok(out)
415    };
416    let low = corners_of(model, bottom)?;
417    let high = corners_of(model, top)?;
418    if low.len() != high.len() {
419        ogeom_bail!(
420            Construction,
421            "lofted sections must have the same corner count, found {} and {}",
422            low.len(),
423            high.len()
424        );
425    }
426    let n = low.len();
427    let centroid = {
428        let mut c = Vector::new(0.0, 0.0, 0.0);
429        for p in low.iter().chain(high.iter()) {
430            c += p.to_vector();
431        }
432        #[allow(clippy::cast_precision_loss)]
433        let count = 2.0 * n as f64;
434        Point::from_vector(c / count)
435    };
436
437    // Shared vertices and edges, then walls and caps referencing them.
438    let vl: Vec<Shape> = low.iter().map(|p| make_vertex(model, *p).shape).collect();
439    let vh: Vec<Shape> = high.iter().map(|p| make_vertex(model, *p).shape).collect();
440    let seg = |model: &mut Model, a: (&Shape, Point), b: (&Shape, Point)| -> OgeomResult<Shape> {
441        let line = LineCurve::segment(a.1, b.1, tol)?;
442        let curve = Curve::Line(line);
443        let domain = curve.domain();
444        Ok(make_edge_between(model, curve, domain, a.0, b.0, tol)?.shape)
445    };
446    let mut low_edges = Vec::with_capacity(n);
447    let mut high_edges = Vec::with_capacity(n);
448    let mut rails = Vec::with_capacity(n);
449    for i in 0..n {
450        let j = (i + 1) % n;
451        low_edges.push(seg(model, (&vl[i], low[i]), (&vl[j], low[j]))?);
452        high_edges.push(seg(model, (&vh[i], high[i]), (&vh[j], high[j]))?);
453        rails.push(seg(model, (&vl[i], low[i]), (&vh[i], high[i]))?);
454    }
455
456    let planar = |model: &mut Model, corners: &[Point], edges: Vec<Shape>| -> OgeomResult<Shape> {
457        let normal = {
458            let mut n = (corners[1] - corners[0]).cross(corners[2] - corners[0]);
459            let m = n.magnitude();
460            if m <= tol.confusion() {
461                ogeom_bail!(Construction, "a loft wall is degenerate");
462            }
463            n /= m;
464            if n.dot(corners[0] - centroid) < 0.0 {
465                -n
466            } else {
467                n
468            }
469        };
470        let skew = corners.iter().any(|p| {
471            Plane::through(
472                corners[0],
473                Direction::new(normal, tol).unwrap_or(Direction::Z),
474            )
475            .distance_to(*p)
476                > tol.confusion() * 10.0
477        });
478        if skew {
479            // A wall between two segments that do not lie in one plane is
480            // the ruled surface between them, and that is exact: bilinear
481            // in the four corners, a B-spline of degree one each way.
482            if corners.len() != 4 || edges.len() != 4 {
483                ogeom_bail!(Construction, "a skew ruled wall has four corners");
484            }
485            return bilinear_wall(model, corners, &edges, centroid, tol);
486        }
487        let plane = Plane::through(corners[0], Direction::new(normal, tol)?);
488        let mut reach = 1.0_f64;
489        for p in corners {
490            reach = reach.max(p.distance(corners[0]) * 2.0);
491        }
492        let surface: SurfaceGeometry =
493            PlaneSurface::over(plane, (-reach, reach), (-reach, reach))?.into();
494        Ok(make_face_with_pcurves(model, surface, &[edges], tol)?.shape)
495    };
496
497    let mut faces: Vec<Shape> = Vec::with_capacity(n + 2);
498    for i in 0..n {
499        let j = (i + 1) % n;
500        faces.push(planar(
501            model,
502            &[low[i], low[j], high[j], high[i]],
503            vec![
504                low_edges[i].clone(),
505                rails[j].clone(),
506                high_edges[i].reversed(),
507                rails[i].reversed(),
508            ],
509        )?);
510    }
511    faces.push(planar(model, &low, low_edges.clone())?);
512    faces.push(planar(model, &high, high_edges.clone())?);
513
514    let sewn = sew(model, &faces, tol)?;
515    if sewn.shells.len() != 1 || !ogeom_algo::is_shell_closed(model, &sewn.shells[0])? {
516        ogeom_bail!(Construction, "the loft did not close");
517    }
518    let mut built = make_solid(model, std::slice::from_ref(&sewn.shells[0]))?;
519    built.history.generate(bottom, built.shape.clone());
520    built.history.generate(top, built.shape.clone());
521    Ok(built)
522}
523
524/// The ruled wall between two straight segments that are not coplanar: the
525/// bilinear patch through its four corners, exact, as a B-spline of degree
526/// one each way.
527///
528/// `corners` run round the wall (low start, low end, high end, high
529/// start) and `edges` walk them in that order, the third and fourth
530/// reversed as the caller's wire has them. Each edge's pcurve is the chart
531/// side it lies on, parameterized by the edge's own range so the two agree
532/// point for point; the wall faces away from `centroid`.
533fn bilinear_wall(
534    model: &mut Model,
535    corners: &[Point],
536    edges: &[Shape],
537    centroid: Point,
538    tol: Tolerances,
539) -> OgeomResult<Shape> {
540    use ogeom_geom::Surface as _;
541    let (p00, p10, p11, p01) = (corners[0], corners[1], corners[2], corners[3]);
542    let grid = ogeom_math::ControlGrid::new(vec![p00, p01, p10, p11], 2, 2)?;
543    let line = ogeom_math::KnotVector::clamped_uniform(1, 2)?;
544    let patch = ogeom_geom::BSplineSurface::new(line.clone(), line, &grid, tol)?;
545    let outward = {
546        let (du, dv) = patch.d1_at(0.5, 0.5, tol)?;
547        let centre = patch.point_at(0.5, 0.5, tol)?;
548        du.cross(dv).dot(centre - centroid) >= 0.0
549    };
550    let surface_id = model
551        .geometry_mut()
552        .add_surface(SurfaceGeometry::BSpline(patch));
553
554    // Each edge's chart side, run in the edge's own direction over the
555    // edge's own parameter range: a degree-one B-spline in the chart,
556    // which is what makes the pcurve the edge's equal parameter for
557    // parameter whatever length the edge has.
558    let sides: [(Point2, Point2); 4] = [
559        (Point2::new(0.0, 0.0), Point2::new(1.0, 0.0)),
560        (Point2::new(1.0, 0.0), Point2::new(1.0, 1.0)),
561        (Point2::new(0.0, 1.0), Point2::new(1.0, 1.0)),
562        (Point2::new(0.0, 0.0), Point2::new(0.0, 1.0)),
563    ];
564    for (edge, (from, to)) in edges.iter().zip(sides) {
565        let range = {
566            let Some(node) = model.node(edge) else {
567                ogeom_bail!(Dangling, "a loft edge is not in this model");
568            };
569            let Some(data) = node.data().as_edge() else {
570                ogeom_bail!(Construction, "a loft edge holds no edge data");
571            };
572            let Some(EdgeRepr::Curve3d { range, .. }) = data.curve3d() else {
573                ogeom_bail!(Construction, "a loft edge has no curve");
574            };
575            *range
576        };
577        let knots = ogeom_math::KnotVector::new(vec![range.0, range.0, range.1, range.1], 1)?;
578        let pcurve = ogeom_geom::BSpline2d::new(knots, vec![from, to], tol)?;
579        ogeom_algo::attach_pcurve(
580            model,
581            edge,
582            pcurve.into(),
583            surface_id,
584            ogeom_topo::Location::identity(),
585            range,
586        )?;
587    }
588    let wire = ogeom_algo::make_wire(model, edges, tol)?.shape;
589    let face = ogeom_algo::make_face_on(model, surface_id, std::slice::from_ref(&wire), tol)?.shape;
590    Ok(if outward { face } else { face.reversed() })
591}
592
593/// A skinned wall and the pieces a caller needs to close it: the rings at
594/// both ends, their exact border curves off the control net, and the chart's
595/// `u` window the ring pcurves span.
596struct SkinnedWall {
597    face: Shape,
598    ring0: Shape,
599    ring1: Shape,
600    curve0: ogeom_geom::Curve,
601    curve1: ogeom_geom::Curve,
602    u_dom: (f64, f64),
603}
604
605/// The wall of a skin over a grid of section samples, closed the way round.
606///
607/// The wall is [`ogeom_geom::fit::fit_surface_grid`]'s surface with each row's
608/// first sample repeated at its end: the row fits pin their ends, so the two
609/// border control columns are *equal* and the seam closes exactly, not
610/// within tolerance. The border iso-curves come straight off the control
611/// net (the v-borders are the fitted sections, planar whenever the
612/// sections are, which is what lets the caps be planes), and every pcurve
613/// is an iso line in the fitted chart, same-parameter by construction.
614fn skinned_wall(
615    model: &mut Model,
616    rows: &[Vec<Point>],
617    shared: (Option<&Shape>, Option<&Shape>),
618    tolerance: f64,
619    tol: Tolerances,
620) -> OgeomResult<SkinnedWall> {
621    use ogeom_geom::Surface as _;
622    let mut closed_rows: Vec<Vec<Point>> = Vec::with_capacity(rows.len());
623    for row in rows {
624        let mut r = row.clone();
625        r.push(row[0]);
626        closed_rows.push(r);
627    }
628    let fitted = ogeom_geom::fit::fit_surface_grid(&closed_rows, 3, tolerance, tol)?;
629    if !fitted.met {
630        ogeom_bail!(
631            NotDone,
632            "the skin reached {} against a target of {tolerance}",
633            fitted.error
634        );
635    }
636    let surface = fitted.curve;
637    let (u_knots, v_knots) = (surface.u_knots().clone(), surface.v_knots().clone());
638    let (k, l, net) = {
639        let grid = surface.grid();
640        let net: Vec<Point> = grid.points().iter().map(|w| (*w).point()).collect();
641        (grid.u_count(), grid.v_count(), net)
642    };
643    let point_at = |i: usize, j: usize| -> Point { net[i * l + j] };
644    let (u_dom, v_dom) = surface.domain();
645
646    // Border curves straight off the net: v-borders are the end sections,
647    // the u-border is the seam.
648    let border_v = |j: usize| -> OgeomResult<ogeom_geom::Curve> {
649        let control: Vec<Point> = (0..k).map(|i| point_at(i, j)).collect();
650        Ok(ogeom_geom::Curve::BSpline(ogeom_geom::BSplineCurve::new(
651            u_knots.clone(),
652            control,
653            tol,
654        )?))
655    };
656    let seam_curve = {
657        let control: Vec<Point> = (0..l).map(|j| point_at(0, j)).collect();
658        ogeom_geom::Curve::BSpline(ogeom_geom::BSplineCurve::new(
659            v_knots.clone(),
660            control,
661            tol,
662        )?)
663    };
664
665    let surface_geo: SurfaceGeometry = surface.into();
666    let surface_id = model.geometry_mut().add_surface(surface_geo.clone());
667
668    // An end ring the neighbouring wall already built is adopted, not
669    // refitted (see `adopt_border`).
670    let slack = fitted.error + tol.confusion();
671    let ring_of = |model: &mut Model,
672                   j: usize,
673                   given: Option<&Shape>|
674     -> OgeomResult<(Shape, ogeom_geom::Curve)> {
675        match given {
676            Some(edge) => {
677                adopt_border(model, edge, &surface_geo, slack, tol)?;
678                Ok((edge.clone(), spine_curve_of(model, edge)?.0))
679            }
680            None => {
681                let curve = border_v(j)?;
682                Ok((make_edge(model, curve.clone(), u_dom, tol)?.shape, curve))
683            }
684        }
685    };
686    let (ring0, curve0) = ring_of(model, 0, shared.0)?;
687    let (ring1, curve1) = ring_of(model, l - 1, shared.1)?;
688    let anchor0 = ogeom_algo::edge_vertices(model, &ring0)?
689        .map(|(a, _)| a)
690        .ok_or_else(|| ogeom_core::ogeom_err!(Construction, "a skinned ring has no vertex"))?;
691    let anchor1 = ogeom_algo::edge_vertices(model, &ring1)?
692        .map(|(a, _)| a)
693        .ok_or_else(|| ogeom_core::ogeom_err!(Construction, "a skinned ring has no vertex"))?;
694    let seam = make_edge_between(model, seam_curve, v_dom, &anchor0, &anchor1, tol)?.shape;
695
696    // Pcurves: rows for the rings, both columns for the seam.
697    let row_line = |v: f64| -> OgeomResult<ogeom_geom::PlanarCurve> {
698        Ok(Line2d::over(
699            ogeom_math::Axis2::new(Point2::new(0.0, v), ogeom_math::Direction2::X),
700            u_dom.0 - 1.0,
701            u_dom.1 + 1.0,
702        )?
703        .into())
704    };
705    let column_line = |u: f64| -> OgeomResult<ogeom_geom::PlanarCurve> {
706        Ok(Line2d::over(
707            ogeom_math::Axis2::new(Point2::new(u, 0.0), ogeom_math::Direction2::Y),
708            v_dom.0 - 1.0,
709            v_dom.1 + 1.0,
710        )?
711        .into())
712    };
713    ogeom_algo::attach_pcurve(
714        model,
715        &ring0,
716        row_line(v_dom.0)?,
717        surface_id,
718        ogeom_topo::Location::identity(),
719        u_dom,
720    )?;
721    ogeom_algo::attach_pcurve(
722        model,
723        &ring1,
724        row_line(v_dom.1)?,
725        surface_id,
726        ogeom_topo::Location::identity(),
727        u_dom,
728    )?;
729    ogeom_algo::attach_seam(
730        model,
731        &seam,
732        column_line(u_dom.0)?,
733        column_line(u_dom.1)?,
734        surface_id,
735        ogeom_topo::Location::identity(),
736        v_dom,
737    )?;
738
739    let wall = {
740        let wire = ogeom_algo::make_wire(
741            model,
742            &[
743                ring0.clone(),
744                seam.clone(),
745                ring1.reversed(),
746                seam.reversed(),
747            ],
748            tol,
749        )?
750        .shape;
751        let face =
752            ogeom_algo::make_face_on(model, surface_id, std::slice::from_ref(&wire), tol)?.shape;
753        // Outward by measurement at the middle of the skin.
754        let mid_u = f64::midpoint(u_dom.0, u_dom.1);
755        let mid_v = f64::midpoint(v_dom.0, v_dom.1);
756        let s_mid = surface_geo.point_at(mid_u, mid_v, tol)?;
757        let (du, dv) = surface_geo.d1_at(mid_u, mid_v, tol)?;
758        let centroid = {
759            let mut c = Vector::new(0.0, 0.0, 0.0);
760            let mut n = 0.0;
761            for row in rows {
762                for p in row {
763                    c += p.to_vector();
764                    n += 1.0;
765                }
766            }
767            Point::from_vector(c / n)
768        };
769        if du.cross(dv).dot(s_mid - centroid) >= 0.0 {
770            face
771        } else {
772            face.reversed()
773        }
774    };
775    Ok(SkinnedWall {
776        face: wall,
777        ring0,
778        ring1,
779        curve0,
780        curve1,
781        u_dom,
782    })
783}
784
785/// A strip closed the *long* way: open across its own width, a smooth loop
786/// along the sweep: one face of a faceted ring, [`skinned_wall`]'s
787/// construction with the chart's roles swapped and the loop made C1 by
788/// [`ogeom_geom::fit::fit_surface_grid_closed_v`]. The rails are the two
789/// closed border loops; the seam is one station's column, used twice.
790fn skinned_ring_strip(
791    model: &mut Model,
792    rows: &[Vec<Point>],
793    outward_hint: Point,
794    shared: [Option<&Shape>; 2],
795    tolerance: f64,
796    tol: Tolerances,
797) -> OgeomResult<(Shape, Shape, Shape)> {
798    use ogeom_geom::Surface as _;
799    // The loop: first row repeated at the end, as the closed fit demands.
800    let mut looped: Vec<Vec<Point>> = rows.to_vec();
801    looped.push(rows[0].clone());
802    let fitted = ogeom_geom::fit::fit_surface_grid_closed_v(&looped, 3, tolerance, tol)?;
803    if !fitted.met {
804        ogeom_bail!(
805            NotDone,
806            "the ring strip reached {} against a target of {tolerance}",
807            fitted.error
808        );
809    }
810    let surface = fitted.curve;
811    let (u_knots, v_knots) = (surface.u_knots().clone(), surface.v_knots().clone());
812    let (k, l, net) = {
813        let grid = surface.grid();
814        let net: Vec<Point> = grid.points().iter().map(|w| (*w).point()).collect();
815        (grid.u_count(), grid.v_count(), net)
816    };
817    let point_at = |i: usize, j: usize| -> Point { net[i * l + j] };
818    let (u_dom, v_dom) = surface.domain();
819
820    // The chart's roles, straight: `u` runs across the strip (open), `v`
821    // around the loop (closed). The rails are v-curves (the closed border
822    // loops at the two u-borders), and the seam is the u-row at the loop's
823    // join, bounding the chart twice as every seam does.
824    let rail_curve = |i: usize| -> OgeomResult<ogeom_geom::Curve> {
825        let control: Vec<Point> = (0..l).map(|j| point_at(i, j)).collect();
826        Ok(ogeom_geom::Curve::BSpline(ogeom_geom::BSplineCurve::new(
827            v_knots.clone(),
828            control,
829            tol,
830        )?))
831    };
832    let seam_curve = {
833        let control: Vec<Point> = (0..k).map(|i| point_at(i, 0)).collect();
834        ogeom_geom::Curve::BSpline(ogeom_geom::BSplineCurve::new(
835            u_knots.clone(),
836            control,
837            tol,
838        )?)
839    };
840    let surface_geo: SurfaceGeometry = surface.into();
841    let surface_id = model.geometry_mut().add_surface(surface_geo.clone());
842
843    // A corner loop shared with the neighbouring strip is one edge for
844    // both: the neighbour built it from its own fit, and this strip's
845    // border is another fit of the same loop, so the edge widens to how
846    // far it honestly sits from this surface. Two independent fits of one
847    // loop can disagree by more than either fit's own error, which is what
848    // the sew refused under a frame that turns fast, and one edge cannot
849    // disagree with itself.
850    let slack = fitted.error + tol.confusion();
851    let rail_of = |model: &mut Model, i: usize, given: Option<&Shape>| -> OgeomResult<Shape> {
852        let Some(edge) = given else {
853            let edge = make_edge(model, rail_curve(i)?, v_dom, tol)?.shape;
854            model.widen(&edge, ogeom_core::Tolerance::new(slack)?)?;
855            return Ok(edge);
856        };
857        let (curve, range) = spine_curve_of(model, edge)?;
858        let mut off: f64 = 0.0;
859        for step in 0..=32 {
860            #[allow(clippy::cast_precision_loss)]
861            let t = range.0 + (range.1 - range.0) * (step as f64) / 32.0;
862            let p = curve.point_at(t, tol)?;
863            off = off.max(ogeom_algo::project_on_surface(&surface_geo, p, 16, tol)?.distance);
864        }
865        model.widen(edge, ogeom_core::Tolerance::new(off + slack)?)?;
866        if let Some((a, b)) = ogeom_algo::edge_vertices(model, edge)? {
867            for v in [&a, &b] {
868                model.widen(v, ogeom_core::Tolerance::new(off + slack)?)?;
869            }
870        }
871        Ok(edge.clone())
872    };
873    let rail0 = rail_of(model, 0, shared[0])?;
874    let rail1 = rail_of(model, k - 1, shared[1])?;
875    let anchor0 = ogeom_algo::edge_vertices(model, &rail0)?
876        .map(|(a, _)| a)
877        .ok_or_else(|| ogeom_core::ogeom_err!(Construction, "a strip rail has no vertex"))?;
878    let anchor1 = ogeom_algo::edge_vertices(model, &rail1)?
879        .map(|(a, _)| a)
880        .ok_or_else(|| ogeom_core::ogeom_err!(Construction, "a strip rail has no vertex"))?;
881    let seam = make_edge_between(model, seam_curve, u_dom, &anchor0, &anchor1, tol)?.shape;
882
883    let row_line = |v: f64| -> OgeomResult<ogeom_geom::PlanarCurve> {
884        Ok(Line2d::over(
885            ogeom_math::Axis2::new(Point2::new(0.0, v), ogeom_math::Direction2::X),
886            u_dom.0 - 1.0,
887            u_dom.1 + 1.0,
888        )?
889        .into())
890    };
891    let column_line = |u: f64| -> OgeomResult<ogeom_geom::PlanarCurve> {
892        Ok(Line2d::over(
893            ogeom_math::Axis2::new(Point2::new(u, 0.0), ogeom_math::Direction2::Y),
894            v_dom.0 - 1.0,
895            v_dom.1 + 1.0,
896        )?
897        .into())
898    };
899    ogeom_algo::attach_pcurve(
900        model,
901        &rail0,
902        column_line(u_dom.0)?,
903        surface_id,
904        ogeom_topo::Location::identity(),
905        v_dom,
906    )?;
907    ogeom_algo::attach_pcurve(
908        model,
909        &rail1,
910        column_line(u_dom.1)?,
911        surface_id,
912        ogeom_topo::Location::identity(),
913        v_dom,
914    )?;
915    ogeom_algo::attach_seam(
916        model,
917        &seam,
918        row_line(v_dom.0)?,
919        row_line(v_dom.1)?,
920        surface_id,
921        ogeom_topo::Location::identity(),
922        u_dom,
923    )?;
924    let wire = ogeom_algo::make_wire(
925        model,
926        &[
927            rail0.clone(),
928            seam.clone(),
929            rail1.reversed(),
930            seam.reversed(),
931        ],
932        tol,
933    )?
934    .shape;
935    let face = ogeom_algo::make_face_on(model, surface_id, std::slice::from_ref(&wire), tol)?.shape;
936    let mid_u = f64::midpoint(u_dom.0, u_dom.1);
937    let mid_v = f64::midpoint(v_dom.0, v_dom.1);
938    let s_mid = surface_geo.point_at(mid_u, mid_v, tol)?;
939    let (du, dv) = surface_geo.d1_at(mid_u, mid_v, tol)?;
940    let face = if du.cross(dv).dot(s_mid - outward_hint) >= 0.0 {
941        face
942    } else {
943        face.reversed()
944    };
945    Ok((face, rail0, rail1))
946}
947
948/// A border edge the neighbouring skin built from its own fit, adopted by
949/// this skin: the edge widens to how far it honestly sits from `surface`
950/// plus this fit's own slack, its vertices with it. Two independent fits of
951/// one row can disagree by more than either fit's own error, and one edge
952/// cannot disagree with itself.
953fn adopt_border(
954    model: &mut Model,
955    edge: &Shape,
956    surface: &SurfaceGeometry,
957    slack: f64,
958    tol: Tolerances,
959) -> OgeomResult<()> {
960    let (curve, range) = spine_curve_of(model, edge)?;
961    let mut off: f64 = 0.0;
962    for step in 0..=32 {
963        #[allow(clippy::cast_precision_loss)]
964        let t = range.0 + (range.1 - range.0) * (step as f64) / 32.0;
965        let p = curve.point_at(t, tol)?;
966        off = off.max(ogeom_algo::project_on_surface(surface, p, 16, tol)?.distance);
967    }
968    model.widen(edge, ogeom_core::Tolerance::new(off + slack)?)?;
969    if let Some((a, b)) = ogeom_algo::edge_vertices(model, edge)? {
970        for v in [&a, &b] {
971            model.widen(v, ogeom_core::Tolerance::new(off + slack)?)?;
972        }
973    }
974    Ok(())
975}
976
977/// An adopted border's image on a surface it was not fitted on: the
978/// border's own points, each read off the surface at its nearest point,
979/// fitted at the border's own parameters so the image is same-parameter
980/// with it. Returned with the parameter range it spans.
981fn adopted_image(
982    model: &Model,
983    edge: &Shape,
984    surface: &SurfaceGeometry,
985    tol: Tolerances,
986) -> OgeomResult<(ogeom_geom::PlanarCurve, (f64, f64))> {
987    const SAMPLES: u32 = 64;
988    let (curve, range) = spine_curve_of(model, edge)?;
989    let mut params = Vec::with_capacity(SAMPLES as usize + 1);
990    let mut image = Vec::with_capacity(SAMPLES as usize + 1);
991    let mut guess: Option<(f64, f64)> = None;
992    for step in 0..=SAMPLES {
993        let t = range.0 + (range.1 - range.0) * f64::from(step) / f64::from(SAMPLES);
994        let p = curve.point_at(t, tol)?;
995        let foot = match guess {
996            Some(g) => ogeom_algo::project_on_surface_from(surface, p, g, tol)?,
997            None => ogeom_algo::project_on_surface(surface, p, 16, tol)?,
998        };
999        guess = Some(foot.parameters);
1000        params.push(t);
1001        image.push(Point2::new(foot.parameters.0, foot.parameters.1));
1002    }
1003    let fitted = ogeom_geom::fit::fit_points_2d_at(&params, &image, 3, tol.confusion(), tol)?;
1004    Ok((fitted.curve.into(), range))
1005}
1006
1007/// A solid skinned over a grid of section samples: [`skinned_wall`] with a
1008/// planar cap over each end ring.
1009/// How a skinned solid's end is closed.
1010#[derive(Debug, Clone, Copy)]
1011enum EndCap {
1012    /// The section is planar: a plane face, its normal pointing out.
1013    Plane(Vector),
1014    /// The section is not: a patch skinned from the ring down to a point
1015    /// inside it, sharing the wall's ring edge.
1016    Skinned,
1017}
1018
1019fn skinned_solid(
1020    model: &mut Model,
1021    rows: &[Vec<Point>],
1022    caps: (EndCap, EndCap),
1023    tolerance: f64,
1024    tol: Tolerances,
1025) -> OgeomResult<Built> {
1026    let wall = skinned_wall(model, rows, (None, None), tolerance, tol)?;
1027    let u_dom = wall.u_dom;
1028    let inside = centroid_of(rows);
1029
1030    let cap = |model: &mut Model,
1031               ring: &Shape,
1032               curve: ogeom_geom::Curve,
1033               outward: Vector|
1034     -> OgeomResult<Shape> {
1035        let at = curve.point_at(u_dom.0, tol)?;
1036        let plane = Plane::through(at, Direction::new(outward, tol)?);
1037        let mut reach = 1.0_f64;
1038        for t in 0..8 {
1039            let p = curve.point_at(u_dom.0 + (u_dom.1 - u_dom.0) * f64::from(t) / 8.0, tol)?;
1040            reach = reach.max(p.distance(at) * 2.0);
1041        }
1042        let cap_surface: SurfaceGeometry =
1043            PlaneSurface::over(plane, (-reach, reach), (-reach, reach))?.into();
1044        let wire = ogeom_algo::make_wire(model, std::slice::from_ref(ring), tol)?.shape;
1045        let face =
1046            ogeom_algo::make_face(model, cap_surface.clone(), std::slice::from_ref(&wire), tol)?
1047                .shape;
1048        let id = {
1049            let Some(node) = model.node(&face) else {
1050                ogeom_bail!(Dangling, "the cap just built is not in this model");
1051            };
1052            let ogeom_topo::NodeData::Face(data) = node.data() else {
1053                ogeom_bail!(Construction, "the cap holds no face data");
1054            };
1055            data.surface
1056        };
1057        let Some(pcurve) = ogeom_intersect::exact_pcurve_of(&curve, &cap_surface, tol) else {
1058            ogeom_bail!(Construction, "a cap edge has no closed-form pcurve");
1059        };
1060        ogeom_algo::attach_pcurve(
1061            model,
1062            ring,
1063            pcurve,
1064            id,
1065            ogeom_topo::Location::identity(),
1066            u_dom,
1067        )?;
1068        Ok(face)
1069    };
1070    let close =
1071        |model: &mut Model, end: EndCap, ring: &Shape, curve: &ogeom_geom::Curve, row: &[Point]| {
1072            match end {
1073                EndCap::Plane(outward) => cap(model, ring, curve.clone(), outward),
1074                EndCap::Skinned => {
1075                    // The ring, a row halfway in, and the point the rest of
1076                    // the ring's rows collapse to: the section's own centroid,
1077                    // which a closed section winds round.
1078                    let apex = centroid_of(std::slice::from_ref(&row.to_vec()));
1079                    let half: Vec<Point> = row
1080                        .iter()
1081                        .map(|p| Point::from_vector((p.to_vector() + apex.to_vector()) * 0.5))
1082                        .collect();
1083                    let rows = [row.to_vec(), half, vec![apex; row.len()]];
1084                    Ok(apex_patch(model, &rows, Some(ring), inside, tolerance, tol)?.0)
1085                }
1086            }
1087        };
1088    let cap0 = close(model, caps.0, &wall.ring0, &wall.curve0, &rows[0])?;
1089    let cap1 = close(
1090        model,
1091        caps.1,
1092        &wall.ring1,
1093        &wall.curve1,
1094        &rows[rows.len() - 1],
1095    )?;
1096
1097    let faces = [wall.face, cap0, cap1];
1098    let sewn = sew(model, &faces, tol)?;
1099    if sewn.shells.len() != 1 || !ogeom_algo::is_shell_closed(model, &sewn.shells[0])? {
1100        ogeom_bail!(Construction, "the skinned solid did not close");
1101    }
1102    make_solid(model, std::slice::from_ref(&sewn.shells[0]))
1103}
1104
1105/// A patch skinned from a ring down to a point: [`skinned_wall`]'s
1106/// construction with the top ring replaced by the apex: a degenerate
1107/// edge on one vertex, bounding the chart's whole top row the way a cone's
1108/// apex bounds a countersink. The ring edge is adopted from `shared`
1109/// where a neighbour already built it, and the face is turned to point
1110/// away from `inside`. Returns the face and its ring edge.
1111fn apex_patch(
1112    model: &mut Model,
1113    rows: &[Vec<Point>],
1114    shared: Option<&Shape>,
1115    inside: Point,
1116    tolerance: f64,
1117    tol: Tolerances,
1118) -> OgeomResult<(Shape, Shape)> {
1119    use ogeom_geom::Surface as _;
1120    let mut closed_rows: Vec<Vec<Point>> = Vec::with_capacity(rows.len());
1121    for row in rows {
1122        let mut r = row.clone();
1123        r.push(row[0]);
1124        closed_rows.push(r);
1125    }
1126    let fitted = ogeom_geom::fit::fit_surface_grid(&closed_rows, 3, tolerance, tol)?;
1127    if !fitted.met {
1128        ogeom_bail!(
1129            NotDone,
1130            "the skin reached {} against a target of {tolerance}",
1131            fitted.error
1132        );
1133    }
1134    let surface = fitted.curve;
1135    let (u_knots, v_knots) = (surface.u_knots().clone(), surface.v_knots().clone());
1136    let (k, l, net) = {
1137        let grid = surface.grid();
1138        let net: Vec<Point> = grid.points().iter().map(|w| (*w).point()).collect();
1139        (grid.u_count(), grid.v_count(), net)
1140    };
1141    let point_at = |i: usize, j: usize| -> Point { net[i * l + j] };
1142    let (u_dom, v_dom) = surface.domain();
1143    let apex = rows[rows.len() - 1][0];
1144
1145    let ring_curve = {
1146        let control: Vec<Point> = (0..k).map(|i| point_at(i, 0)).collect();
1147        ogeom_geom::Curve::BSpline(ogeom_geom::BSplineCurve::new(
1148            u_knots.clone(),
1149            control,
1150            tol,
1151        )?)
1152    };
1153    let seam_curve = {
1154        let control: Vec<Point> = (0..l).map(|j| point_at(0, j)).collect();
1155        ogeom_geom::Curve::BSpline(ogeom_geom::BSplineCurve::new(
1156            v_knots.clone(),
1157            control,
1158            tol,
1159        )?)
1160    };
1161    let surface_geo: SurfaceGeometry = surface.into();
1162    let surface_id = model.geometry_mut().add_surface(surface_geo.clone());
1163
1164    // The ring a neighbour built is adopted, not refitted (see `adopt_border`).
1165    let ring0 = match shared {
1166        Some(edge) => {
1167            adopt_border(
1168                model,
1169                edge,
1170                &surface_geo,
1171                fitted.error + tol.confusion(),
1172                tol,
1173            )?;
1174            edge.clone()
1175        }
1176        None => make_edge(model, ring_curve.clone(), u_dom, tol)?.shape,
1177    };
1178    let anchor0 = ogeom_algo::edge_vertices(model, &ring0)?
1179        .map(|(a, _)| a)
1180        .ok_or_else(|| ogeom_core::ogeom_err!(Construction, "a skinned ring has no vertex"))?;
1181    let apex_vertex = model.add_vertex(VertexData::new(apex));
1182    let apex_edge = {
1183        let mut data = EdgeData::new();
1184        data.degenerate = true;
1185        model.add_edge(data, &[apex_vertex.clone(), apex_vertex.clone()])?
1186    };
1187    let seam = make_edge_between(model, seam_curve, v_dom, &anchor0, &apex_vertex, tol)?.shape;
1188
1189    let row_line = |v: f64| -> OgeomResult<ogeom_geom::PlanarCurve> {
1190        Ok(Line2d::over(
1191            ogeom_math::Axis2::new(Point2::new(0.0, v), ogeom_math::Direction2::X),
1192            u_dom.0 - 1.0,
1193            u_dom.1 + 1.0,
1194        )?
1195        .into())
1196    };
1197    let column_line = |u: f64| -> OgeomResult<ogeom_geom::PlanarCurve> {
1198        Ok(Line2d::over(
1199            ogeom_math::Axis2::new(Point2::new(u, 0.0), ogeom_math::Direction2::Y),
1200            v_dom.0 - 1.0,
1201            v_dom.1 + 1.0,
1202        )?
1203        .into())
1204    };
1205    ogeom_algo::attach_pcurve(
1206        model,
1207        &ring0,
1208        row_line(v_dom.0)?,
1209        surface_id,
1210        ogeom_topo::Location::identity(),
1211        u_dom,
1212    )?;
1213    // The apex bounds the chart's whole top row while covering no distance:
1214    // the degenerate edge carries the row's pcurve, exactly as a cone's apex
1215    // does after the reader synthesises it.
1216    ogeom_algo::attach_pcurve(
1217        model,
1218        &apex_edge,
1219        row_line(v_dom.1)?,
1220        surface_id,
1221        ogeom_topo::Location::identity(),
1222        u_dom,
1223    )?;
1224    ogeom_algo::attach_seam(
1225        model,
1226        &seam,
1227        column_line(u_dom.0)?,
1228        column_line(u_dom.1)?,
1229        surface_id,
1230        ogeom_topo::Location::identity(),
1231        v_dom,
1232    )?;
1233
1234    let wire = ogeom_algo::make_wire(
1235        model,
1236        &[
1237            ring0.clone(),
1238            seam.clone(),
1239            apex_edge.reversed(),
1240            seam.reversed(),
1241        ],
1242        tol,
1243    )?
1244    .shape;
1245    let face = ogeom_algo::make_face_on(model, surface_id, std::slice::from_ref(&wire), tol)?.shape;
1246    let mid_u = f64::midpoint(u_dom.0, u_dom.1);
1247    let mid_v = f64::midpoint(v_dom.0, v_dom.1);
1248    let s_mid = surface_geo.point_at(mid_u, mid_v, tol)?;
1249    let (du, dv) = surface_geo.d1_at(mid_u, mid_v, tol)?;
1250    let face = if du.cross(dv).dot(s_mid - inside) >= 0.0 {
1251        face
1252    } else {
1253        face.reversed()
1254    };
1255    Ok((face, ring0))
1256}
1257
1258/// The mean of every point in every row.
1259fn centroid_of(rows: &[Vec<Point>]) -> Point {
1260    let mut c = Vector::new(0.0, 0.0, 0.0);
1261    let mut n = 0.0;
1262    for row in rows {
1263        for p in row {
1264            c += p.to_vector();
1265            n += 1.0;
1266        }
1267    }
1268    Point::from_vector(c / n)
1269}
1270
1271/// A solid skinned down to a point: [`apex_patch`] for the wall, and one
1272/// cap at the open end; the apex end closes by construction.
1273fn skinned_solid_to_apex(
1274    model: &mut Model,
1275    rows: &[Vec<Point>],
1276    cap_outward: Vector,
1277    tolerance: f64,
1278    tol: Tolerances,
1279) -> OgeomResult<Built> {
1280    use ogeom_geom::Curve3d as _;
1281    let inside = centroid_of(rows);
1282    let (wall, ring0) = apex_patch(model, rows, None, inside, tolerance, tol)?;
1283    let (ring_curve, u_dom) = {
1284        let (curve, range) = spine_curve_of(model, &ring0)?;
1285        (curve, range)
1286    };
1287
1288    // One cap, on the open end; the machinery is skinned_solid's, inlined
1289    // for the single ring.
1290    let cap = {
1291        let at = ring_curve.point_at(u_dom.0, tol)?;
1292        let plane = Plane::through(at, Direction::new(cap_outward, tol)?);
1293        let mut reach = 1.0_f64;
1294        for t in 0..8 {
1295            let p = ring_curve.point_at(u_dom.0 + (u_dom.1 - u_dom.0) * f64::from(t) / 8.0, tol)?;
1296            reach = reach.max(p.distance(at) * 2.0);
1297        }
1298        let cap_surface: SurfaceGeometry =
1299            PlaneSurface::over(plane, (-reach, reach), (-reach, reach))?.into();
1300        let wire = ogeom_algo::make_wire(model, std::slice::from_ref(&ring0), tol)?.shape;
1301        let face =
1302            ogeom_algo::make_face(model, cap_surface.clone(), std::slice::from_ref(&wire), tol)?
1303                .shape;
1304        let id = {
1305            let Some(node) = model.node(&face) else {
1306                ogeom_bail!(Dangling, "the cap just built is not in this model");
1307            };
1308            let ogeom_topo::NodeData::Face(data) = node.data() else {
1309                ogeom_bail!(Construction, "the cap holds no face data");
1310            };
1311            data.surface
1312        };
1313        let Some(pcurve) = ogeom_intersect::exact_pcurve_of(&ring_curve, &cap_surface, tol) else {
1314            ogeom_bail!(Construction, "a cap edge has no closed-form pcurve");
1315        };
1316        ogeom_algo::attach_pcurve(
1317            model,
1318            &ring0,
1319            pcurve,
1320            id,
1321            ogeom_topo::Location::identity(),
1322            u_dom,
1323        )?;
1324        face
1325    };
1326
1327    let faces = [wall, cap];
1328    let sewn = sew(model, &faces, tol)?;
1329    if sewn.shells.len() != 1 || !ogeom_algo::is_shell_closed(model, &sewn.shells[0])? {
1330        ogeom_bail!(Construction, "the skinned apex solid did not close");
1331    }
1332    make_solid(model, std::slice::from_ref(&sewn.shells[0]))
1333}
1334
1335/// A solid skinned over a grid of sections that loops back on itself: the
1336/// wall is one face closed in both chart directions, no caps at all.
1337///
1338/// The `u` seam closes the way every skin's does (pinned row ends), and
1339/// the `v` loop closes through [`ogeom_geom::fit::fit_surface_grid_closed_v`],
1340/// C1 across the join. All four boundary traversals are two seam edges used
1341/// twice, anchored at one shared vertex, exactly as a torus bounds itself.
1342fn closed_skinned_solid(
1343    model: &mut Model,
1344    rows: &[Vec<Point>],
1345    tolerance: f64,
1346    tol: Tolerances,
1347) -> OgeomResult<Built> {
1348    let shell = closed_skinned_shell(model, rows, tolerance, tol)?;
1349    make_solid(model, std::slice::from_ref(&shell))
1350}
1351
1352/// The closed skin as a shell, for callers assembling solids with voids;
1353/// a holed profile's ring is one outer shell and one per tunnel.
1354fn closed_skinned_shell(
1355    model: &mut Model,
1356    rows: &[Vec<Point>],
1357    tolerance: f64,
1358    tol: Tolerances,
1359) -> OgeomResult<Shape> {
1360    use ogeom_geom::Surface as _;
1361    let mut closed_rows: Vec<Vec<Point>> = Vec::with_capacity(rows.len() + 1);
1362    for row in rows {
1363        let mut r = row.clone();
1364        r.push(row[0]);
1365        closed_rows.push(r);
1366    }
1367    closed_rows.push(closed_rows[0].clone());
1368    let fitted = ogeom_geom::fit::fit_surface_grid_closed_v(&closed_rows, 3, tolerance, tol)?;
1369    if !fitted.met {
1370        ogeom_bail!(
1371            NotDone,
1372            "the closed skin reached {} against a target of {tolerance}",
1373            fitted.error
1374        );
1375    }
1376    let surface = fitted.curve;
1377    let (u_knots, v_knots) = (surface.u_knots().clone(), surface.v_knots().clone());
1378    let (k, l, net) = {
1379        let grid = surface.grid();
1380        let net: Vec<Point> = grid.points().iter().map(|w| (*w).point()).collect();
1381        (grid.u_count(), grid.v_count(), net)
1382    };
1383    let point_at = |i: usize, j: usize| -> Point { net[i * l + j] };
1384    let (u_dom, v_dom) = surface.domain();
1385
1386    // Both seams straight off the net: the u-run at v's join, and the v-run
1387    // at u's.
1388    let along_u = {
1389        let control: Vec<Point> = (0..k).map(|i| point_at(i, 0)).collect();
1390        ogeom_geom::Curve::BSpline(ogeom_geom::BSplineCurve::new(u_knots, control, tol)?)
1391    };
1392    let along_v = {
1393        let control: Vec<Point> = (0..l).map(|j| point_at(0, j)).collect();
1394        ogeom_geom::Curve::BSpline(ogeom_geom::BSplineCurve::new(v_knots, control, tol)?)
1395    };
1396    let surface_geo: SurfaceGeometry = surface.into();
1397    let surface_id = model.geometry_mut().add_surface(surface_geo.clone());
1398
1399    let u_edge = make_edge(model, along_u, u_dom, tol)?.shape;
1400    let Some((corner, _)) = ogeom_algo::edge_vertices(model, &u_edge)? else {
1401        ogeom_bail!(Construction, "the closed skin's seam has no vertex");
1402    };
1403    let v_edge = make_edge_between(model, along_v, v_dom, &corner, &corner, tol)?.shape;
1404
1405    let row_line = |v: f64| -> OgeomResult<ogeom_geom::PlanarCurve> {
1406        Ok(Line2d::over(
1407            ogeom_math::Axis2::new(Point2::new(0.0, v), ogeom_math::Direction2::X),
1408            u_dom.0 - 1.0,
1409            u_dom.1 + 1.0,
1410        )?
1411        .into())
1412    };
1413    let column_line = |u: f64| -> OgeomResult<ogeom_geom::PlanarCurve> {
1414        Ok(Line2d::over(
1415            ogeom_math::Axis2::new(Point2::new(u, 0.0), ogeom_math::Direction2::Y),
1416            v_dom.0 - 1.0,
1417            v_dom.1 + 1.0,
1418        )?
1419        .into())
1420    };
1421    // The u-run is a seam in v (the same curve at both rows), and the
1422    // v-run a seam in u.
1423    ogeom_algo::attach_seam(
1424        model,
1425        &u_edge,
1426        row_line(v_dom.0)?,
1427        row_line(v_dom.1)?,
1428        surface_id,
1429        ogeom_topo::Location::identity(),
1430        u_dom,
1431    )?;
1432    ogeom_algo::attach_seam(
1433        model,
1434        &v_edge,
1435        column_line(u_dom.1)?,
1436        column_line(u_dom.0)?,
1437        surface_id,
1438        ogeom_topo::Location::identity(),
1439        v_dom,
1440    )?;
1441
1442    let wire = ogeom_algo::make_wire(
1443        model,
1444        &[
1445            u_edge.clone(),
1446            v_edge.clone(),
1447            u_edge.reversed(),
1448            v_edge.reversed(),
1449        ],
1450        tol,
1451    )?
1452    .shape;
1453    let face = ogeom_algo::make_face_on(model, surface_id, std::slice::from_ref(&wire), tol)?.shape;
1454    let centroid = {
1455        let mut c = Vector::new(0.0, 0.0, 0.0);
1456        let mut n = 0.0;
1457        for row in rows {
1458            for p in row {
1459                c += p.to_vector();
1460                n += 1.0;
1461            }
1462        }
1463        Point::from_vector(c / n)
1464    };
1465    let mid_u = f64::midpoint(u_dom.0, u_dom.1);
1466    let mid_v = f64::midpoint(v_dom.0, v_dom.1);
1467    let s_mid = surface_geo.point_at(mid_u, mid_v, tol)?;
1468    let (du, dv) = surface_geo.d1_at(mid_u, mid_v, tol)?;
1469    let face = if du.cross(dv).dot(s_mid - centroid) >= 0.0 {
1470        face
1471    } else {
1472        face.reversed()
1473    };
1474
1475    let sewn = sew(model, std::slice::from_ref(&face), tol)?;
1476    if sewn.shells.len() != 1 || !ogeom_algo::is_shell_closed(model, &sewn.shells[0])? {
1477        ogeom_bail!(Construction, "the closed skin did not close");
1478    }
1479    Ok(sewn.shells[0].clone())
1480}
1481
1482/// The loft to a point: a wire section closing onto a single apex vertex.
1483///
1484/// A circle takes the cone the revolved primitives already build, apex on
1485/// its axis or refused; a polygon takes exact planar triangle walls, sound
1486/// for *any* apex: a skew pyramid's walls are still triangles.
1487fn loft_to_point(
1488    model: &mut Model,
1489    section: &Shape,
1490    apex: &Shape,
1491    tol: Tolerances,
1492) -> OgeomResult<Built> {
1493    if !ogeom_algo::is_wire_closed(model, section, tol)? {
1494        ogeom_bail!(Construction, "a loft section must be closed");
1495    }
1496    let apex_point = {
1497        let Some(data) = model.node(apex).and_then(|n| n.data().as_vertex()) else {
1498            ogeom_bail!(Construction, "the apex vertex holds no data");
1499        };
1500        data.point
1501    };
1502
1503    // The circular case: a cone, apex on the axis.
1504    let edges = explore(model, section, Filter::OfType(ShapeType::Edge))?;
1505    if edges.len() == 1
1506        && let Some(data) = model.node(&edges[0]).and_then(|n| n.data().as_edge())
1507        && let Some(EdgeRepr::Curve3d { curve, .. }) = data.curve3d()
1508        && let Some(Curve::Circle(c)) = model.geometry().curve(*curve)
1509    {
1510        let circle = c.circle();
1511        let axis = circle.frame().z().vector();
1512        let rise = apex_point - circle.centre();
1513        let height = rise.dot(axis);
1514        if rise.cross(axis).magnitude() > tol.confusion() * 10.0 {
1515            ogeom_bail!(
1516                Construction,
1517                "a circle lofts to a point on its own axis; the oblique cone \
1518                 needs the skinned machinery; see docs/PARITY.md, offset.loft"
1519            );
1520        }
1521        if height.abs() <= tol.confusion() {
1522            ogeom_bail!(Construction, "the apex sits in the section's own plane");
1523        }
1524        let base = if height > 0.0 {
1525            circle.frame()
1526        } else {
1527            Frame::new(
1528                circle.centre(),
1529                -circle.frame().z(),
1530                circle.frame().x(),
1531                tol,
1532            )?
1533        };
1534        let mut built =
1535            ogeom_algo::make_cone(model, base, circle.radius(), 0.0, height.abs(), tol)?;
1536        built.history.generate(section, built.shape.clone());
1537        built.history.generate(apex, built.shape.clone());
1538        return Ok(built);
1539    }
1540
1541    // The polygonal case: exact triangle walls to a shared apex.
1542    let mut corners: Vec<Point> = Vec::new();
1543    for edge in model.ordered_children_of(section)? {
1544        let Some(data) = model.node(&edge).and_then(|n| n.data().as_edge()) else {
1545            ogeom_bail!(Construction, "a section edge holds no data");
1546        };
1547        let Some(EdgeRepr::Curve3d { curve, range, .. }) = data.curve3d() else {
1548            ogeom_bail!(Construction, "a section edge has no curve");
1549        };
1550        let Some(Curve::Line(line)) = model.geometry().curve(*curve).cloned() else {
1551            ogeom_bail!(
1552                Construction,
1553                "a mixed or curved section lofts to a point through the \
1554                 skinned machinery; see docs/PARITY.md, offset.loft"
1555            );
1556        };
1557        let t = if edge.orientation() == ogeom_topo::Orientation::Reversed {
1558            range.1
1559        } else {
1560            range.0
1561        };
1562        corners.push(ogeom_geom::Curve::Line(line).point_at(t, tol)?);
1563    }
1564    if corners.len() < 3 {
1565        ogeom_bail!(Construction, "a pyramid needs at least three base corners");
1566    }
1567    let apex_vertex = ogeom_algo::make_vertex(model, apex_point).shape;
1568    let base_vertices: Vec<Shape> = corners
1569        .iter()
1570        .map(|p| ogeom_algo::make_vertex(model, *p).shape)
1571        .collect();
1572    let segment =
1573        |model: &mut Model, from: (&Shape, Point), to: (&Shape, Point)| -> OgeomResult<Shape> {
1574            let line = ogeom_geom::LineCurve::segment(from.1, to.1, tol)?;
1575            let curve: Curve = line.into();
1576            let domain = curve.domain();
1577            Ok(make_edge_between(model, curve, domain, from.0, to.0, tol)?.shape)
1578        };
1579    let count = corners.len();
1580    let mut base_edges = Vec::with_capacity(count);
1581    let mut rails = Vec::with_capacity(count);
1582    for i in 0..count {
1583        let next = (i + 1) % count;
1584        base_edges.push(segment(
1585            model,
1586            (&base_vertices[i], corners[i]),
1587            (&base_vertices[next], corners[next]),
1588        )?);
1589        rails.push(segment(
1590            model,
1591            (&base_vertices[i], corners[i]),
1592            (&apex_vertex, apex_point),
1593        )?);
1594    }
1595    let centroid = {
1596        let mut c = Vector::new(0.0, 0.0, 0.0);
1597        for p in &corners {
1598            c += p.to_vector();
1599        }
1600        #[allow(clippy::cast_precision_loss)]
1601        Point::from_vector(c / count as f64 / 4.0 * 3.0 + apex_point.to_vector() / 4.0)
1602    };
1603    let planar = |model: &mut Model, pts: [Point; 3], walk: Vec<Shape>| -> OgeomResult<Shape> {
1604        let n = (pts[1] - pts[0]).cross(pts[2] - pts[0]);
1605        let m = n.magnitude();
1606        if m <= tol.confusion() {
1607            ogeom_bail!(Construction, "a wall of the pyramid is degenerate");
1608        }
1609        let mut outward = n / m;
1610        if outward.dot(pts[0] - centroid) < 0.0 {
1611            outward = -outward;
1612        }
1613        let plane = ogeom_math::Plane::through(pts[0], Direction::new(outward, tol)?);
1614        let mut reach = 1.0_f64;
1615        for p in pts {
1616            reach = reach.max(p.distance(pts[0]) * 2.0);
1617        }
1618        let surface: SurfaceGeometry =
1619            PlaneSurface::over(plane, (-reach, reach), (-reach, reach))?.into();
1620        let id = model.geometry_mut().add_surface(surface.clone());
1621        let signed = {
1622            let (du, dv) = {
1623                use ogeom_geom::Surface as _;
1624                surface.d1_at(0.0, 0.0, tol)?
1625            };
1626            du.cross(dv).dot(outward) >= 0.0
1627        };
1628        let mut wired = Vec::with_capacity(walk.len());
1629        for used in &walk {
1630            let (curve, range) = spine_curve_of(model, used)?;
1631            let Some(pcurve) = ogeom_intersect::exact_pcurve_of(&curve, &surface, tol) else {
1632                ogeom_bail!(Construction, "a wall edge has no closed-form pcurve");
1633            };
1634            ogeom_algo::attach_pcurve(
1635                model,
1636                used,
1637                pcurve,
1638                id,
1639                ogeom_topo::Location::identity(),
1640                range,
1641            )?;
1642            wired.push(used.clone());
1643        }
1644        let wire = ogeom_algo::make_wire(model, &wired, tol)?.shape;
1645        let face = ogeom_algo::make_face_on(model, id, std::slice::from_ref(&wire), tol)?.shape;
1646        Ok(if signed { face } else { face.reversed() })
1647    };
1648    let mut faces = Vec::with_capacity(count + 1);
1649    for i in 0..count {
1650        let next = (i + 1) % count;
1651        faces.push(planar(
1652            model,
1653            [corners[i], corners[next], apex_point],
1654            vec![
1655                base_edges[i].clone(),
1656                rails[next].clone(),
1657                rails[i].reversed(),
1658            ],
1659        )?);
1660    }
1661    // The base cap: all corners, wound against the walls.
1662    let base_walk: Vec<Shape> = (0..count).rev().map(|i| base_edges[i].reversed()).collect();
1663    faces.push({
1664        let n = (corners[1] - corners[0]).cross(corners[2] - corners[0]);
1665        let mut outward = n / n.magnitude();
1666        if outward.dot(corners[0] - centroid) < 0.0 {
1667            outward = -outward;
1668        }
1669        let plane = ogeom_math::Plane::through(corners[0], Direction::new(outward, tol)?);
1670        let mut reach = 1.0_f64;
1671        for p in &corners {
1672            reach = reach.max(p.distance(corners[0]) * 2.0);
1673        }
1674        let surface: SurfaceGeometry =
1675            PlaneSurface::over(plane, (-reach, reach), (-reach, reach))?.into();
1676        let id = model.geometry_mut().add_surface(surface.clone());
1677        for used in &base_walk {
1678            let (curve, range) = spine_curve_of(model, used)?;
1679            let Some(pcurve) = ogeom_intersect::exact_pcurve_of(&curve, &surface, tol) else {
1680                ogeom_bail!(Construction, "a base edge has no closed-form pcurve");
1681            };
1682            ogeom_algo::attach_pcurve(
1683                model,
1684                used,
1685                pcurve,
1686                id,
1687                ogeom_topo::Location::identity(),
1688                range,
1689            )?;
1690        }
1691        let wire = ogeom_algo::make_wire(model, &base_walk, tol)?.shape;
1692        let face = ogeom_algo::make_face_on(model, id, std::slice::from_ref(&wire), tol)?.shape;
1693        let signed = {
1694            use ogeom_geom::Surface as _;
1695            let (du, dv) = surface.d1_at(0.0, 0.0, tol)?;
1696            du.cross(dv).dot(outward) >= 0.0
1697        };
1698        if signed { face } else { face.reversed() }
1699    });
1700
1701    let sewn = sew(model, &faces, tol)?;
1702    if sewn.shells.len() != 1 || !ogeom_algo::is_shell_closed(model, &sewn.shells[0])? {
1703        ogeom_bail!(Construction, "the pyramid did not close");
1704    }
1705    let mut built = make_solid(model, std::slice::from_ref(&sewn.shells[0]))?;
1706    built.history.generate(section, built.shape.clone());
1707    built.history.generate(apex, built.shape.clone());
1708    Ok(built)
1709}
1710
1711/// Loft through sections with the start of each row named by the caller.
1712///
1713/// [`make_loft_skinned`] leaves alignment to each section's own traversal
1714/// start; this sibling takes one hint per section (a point near where its
1715/// row should begin) and rotates each sampling there, which is how a
1716/// caller untwists a loft whose wires happen to start in different places.
1717///
1718/// # Errors
1719///
1720/// As [`make_loft_skinned`], and additionally if the hints do not pair up
1721/// with the sections.
1722pub fn make_loft_skinned_aligned(
1723    model: &mut Model,
1724    sections: &[Shape],
1725    hints: &[Point],
1726    tolerance: f64,
1727    tol: Tolerances,
1728) -> OgeomResult<Built> {
1729    if hints.len() != sections.len() {
1730        ogeom_bail!(
1731            Construction,
1732            "{} hints against {} sections; each section names its own start",
1733            hints.len(),
1734            sections.len()
1735        );
1736    }
1737    if sections.len() < 2 {
1738        ogeom_bail!(Construction, "a loft needs at least two sections");
1739    }
1740    const AROUND: usize = 48;
1741    let mut rows: Vec<Vec<Point>> = Vec::with_capacity(sections.len());
1742    let mut planes: Vec<Plane> = Vec::with_capacity(sections.len());
1743    for (wire, hint) in sections.iter().zip(hints) {
1744        if model.kind_of(wire)? != ShapeType::Wire {
1745            ogeom_bail!(Construction, "a loft section is a closed wire");
1746        }
1747        if !ogeom_algo::is_wire_closed(model, wire, tol)? {
1748            ogeom_bail!(Construction, "a loft section must be closed");
1749        }
1750        let Some(plane) = ogeom_algo::find_plane(model, wire, tol)? else {
1751            ogeom_bail!(Construction, "a loft section must be planar");
1752        };
1753        planes.push(plane);
1754        rows.push(sample_wire_from(model, wire, AROUND, Some(*hint), tol)?);
1755    }
1756    let outward0 = {
1757        let towards = rows[1][0] - rows[0][0];
1758        let n = planes[0].normal().vector();
1759        if n.dot(towards) > 0.0 { -n } else { n }
1760    };
1761    let outward1 = {
1762        let towards = rows[rows.len() - 2][0] - rows[rows.len() - 1][0];
1763        let n = planes[planes.len() - 1].normal().vector();
1764        if n.dot(towards) > 0.0 { -n } else { n }
1765    };
1766    let mut built = skinned_solid(
1767        model,
1768        &rows,
1769        (EndCap::Plane(outward0), EndCap::Plane(outward1)),
1770        tolerance,
1771        tol,
1772    )?;
1773    for section in sections {
1774        built.history.generate(section, built.shape.clone());
1775    }
1776    Ok(built)
1777}
1778
1779/// Loft a ring through closed planar sections that loop back to the first.
1780///
1781/// [`make_loft_skinned`]'s closed sibling: the sections are sampled the same
1782/// way, the skin runs through all of them and back to the start, C1 across
1783/// the loop, and there are no caps: the result bounds itself the way a
1784/// torus does. The sections are *not* repeated: the loop-back is the
1785/// construction's own.
1786///
1787/// The closed join costs freedom: a sparse loop fits only loosely, and the
1788/// refusal quotes the deviation it honestly reached. A loop that wants a
1789/// tight tolerance wants sections dense enough to bend around: in
1790/// practice, a dozen and up.
1791///
1792/// # Errors
1793///
1794/// As [`make_loft_skinned`], needing at least three sections;
1795/// [`OgeomError::NotDone`](ogeom_core::OgeomError::NotDone) if the closed
1796/// skin cannot reach the tolerance.
1797pub fn make_loft_skinned_closed(
1798    model: &mut Model,
1799    sections: &[Shape],
1800    tolerance: f64,
1801    tol: Tolerances,
1802) -> OgeomResult<Built> {
1803    if sections.len() < 3 {
1804        ogeom_bail!(Construction, "a closed loft needs at least three sections");
1805    }
1806    const AROUND: usize = 48;
1807    let mut rows: Vec<Vec<Point>> = Vec::with_capacity(sections.len());
1808    for wire in sections {
1809        if model.kind_of(wire)? != ShapeType::Wire {
1810            ogeom_bail!(Construction, "a loft section is a closed wire");
1811        }
1812        if !ogeom_algo::is_wire_closed(model, wire, tol)? {
1813            ogeom_bail!(Construction, "a loft section must be closed");
1814        }
1815        if ogeom_algo::find_plane(model, wire, tol)?.is_none() {
1816            ogeom_bail!(Construction, "a loft section must be planar");
1817        }
1818        rows.push(sample_wire(model, wire, AROUND, tol)?);
1819    }
1820    let mut built = closed_skinned_solid(model, &rows, tolerance, tol)?;
1821    for section in sections {
1822        built.history.generate(section, built.shape.clone());
1823    }
1824    Ok(built)
1825}
1826
1827/// A skinned strip: one open patch of a sweep, with its border edges.
1828///
1829/// The wall of a *faceted* profile cannot be one closed skin (a fit cannot
1830/// speak a corner), so each profile edge sweeps its own strip, cornered at
1831/// the caller's shared vertices, and the strips weld along their rails by
1832/// the tolerance the fit honestly carries.
1833struct SkinnedStrip {
1834    face: Shape,
1835    /// The border along the first station, from `corners.0` to `corners.1`.
1836    bottom: Shape,
1837    /// The border along the last station, from `corners.2` to `corners.3`.
1838    top: Shape,
1839    /// The rail along the profile edge's start, from `corners.0` to `corners.2`.
1840    rail0: Shape,
1841    /// The rail along the profile edge's end, from `corners.1` to `corners.3`.
1842    rail1: Shape,
1843}
1844
1845/// Skin an open grid of samples (stations by profile-edge samples) into
1846/// one strip. `corners` are the caller's vertices at (first station, edge
1847/// start), (first, end), (last, start), (last, end), shared with the
1848/// neighbouring strips so the wires chain. `shared` are borders a
1849/// neighbouring strip already built (bottom, top, start rail, end rail),
1850/// adopted as they are (see `adopt_border`).
1851#[allow(clippy::too_many_arguments, reason = "one strip, spelled out")]
1852fn skinned_strip(
1853    model: &mut Model,
1854    rows: &[Vec<Point>],
1855    corners: (&Shape, &Shape, &Shape, &Shape),
1856    shared: [Option<&Shape>; 4],
1857    outward_hint: Point,
1858    hole: bool,
1859    tolerance: f64,
1860    tol: Tolerances,
1861) -> OgeomResult<SkinnedStrip> {
1862    use ogeom_geom::Surface as _;
1863    // A strip whose every row lies in one plane (a straight profile edge
1864    // down a straight run, a flat face of the profile along a planar
1865    // spine) is that plane, exactly: a coplanar neighbour then melts with
1866    // it on the one surface two fits of it would never agree on.
1867    if let Some(plane) = plane_of_rows(rows, tol) {
1868        return planar_strip(
1869            model,
1870            rows,
1871            plane,
1872            corners,
1873            shared,
1874            outward_hint,
1875            hole,
1876            tolerance,
1877            tol,
1878        );
1879    }
1880    let fitted = ogeom_geom::fit::fit_surface_grid(rows, 3, tolerance, tol)?;
1881    if !fitted.met {
1882        ogeom_bail!(
1883            NotDone,
1884            "the strip reached {} against a target of {tolerance}",
1885            fitted.error
1886        );
1887    }
1888    let error = fitted.error.max(tol.confusion());
1889    let surface = fitted.curve;
1890    let (u_knots, v_knots) = (surface.u_knots().clone(), surface.v_knots().clone());
1891    let (k, l, net) = {
1892        let grid = surface.grid();
1893        let net: Vec<Point> = grid.points().iter().map(|w| (*w).point()).collect();
1894        (grid.u_count(), grid.v_count(), net)
1895    };
1896    let point_at = |i: usize, j: usize| -> Point { net[i * l + j] };
1897    let (u_dom, v_dom) = surface.domain();
1898
1899    let u_curve = |j: usize| -> OgeomResult<ogeom_geom::Curve> {
1900        let control: Vec<Point> = (0..k).map(|i| point_at(i, j)).collect();
1901        Ok(ogeom_geom::Curve::BSpline(ogeom_geom::BSplineCurve::new(
1902            u_knots.clone(),
1903            control,
1904            tol,
1905        )?))
1906    };
1907    let v_curve = |i: usize| -> OgeomResult<ogeom_geom::Curve> {
1908        let control: Vec<Point> = (0..l).map(|j| point_at(i, j)).collect();
1909        Ok(ogeom_geom::Curve::BSpline(ogeom_geom::BSplineCurve::new(
1910            v_knots.clone(),
1911            control,
1912            tol,
1913        )?))
1914    };
1915    let surface_geo: SurfaceGeometry = surface.into();
1916    let surface_id = model.geometry_mut().add_surface(surface_geo.clone());
1917
1918    let (c00, c10, c01, c11) = corners;
1919    // A border the neighbouring strip already built is adopted, not
1920    // refitted (see `adopt_border`).
1921    let border = |model: &mut Model,
1922                  given: Option<&Shape>,
1923                  curve: ogeom_geom::Curve,
1924                  range: (f64, f64),
1925                  from: &Shape,
1926                  to: &Shape|
1927     -> OgeomResult<Shape> {
1928        match given {
1929            Some(edge) => {
1930                adopt_border(model, edge, &surface_geo, error + tol.confusion(), tol)?;
1931                Ok(edge.clone())
1932            }
1933            None => Ok(make_edge_between(model, curve, range, from, to, tol)?.shape),
1934        }
1935    };
1936    let bottom = border(model, shared[0], u_curve(0)?, u_dom, c00, c10)?;
1937    let top = border(model, shared[1], u_curve(l - 1)?, u_dom, c01, c11)?;
1938    let rail0 = border(model, shared[2], v_curve(0)?, v_dom, c00, c01)?;
1939    let rail1 = border(model, shared[3], v_curve(k - 1)?, v_dom, c10, c11)?;
1940
1941    let row_line = |v: f64| -> OgeomResult<ogeom_geom::PlanarCurve> {
1942        Ok(Line2d::over(
1943            ogeom_math::Axis2::new(Point2::new(0.0, v), ogeom_math::Direction2::X),
1944            u_dom.0 - 1.0,
1945            u_dom.1 + 1.0,
1946        )?
1947        .into())
1948    };
1949    let column_line = |u: f64| -> OgeomResult<ogeom_geom::PlanarCurve> {
1950        Ok(Line2d::over(
1951            ogeom_math::Axis2::new(Point2::new(u, 0.0), ogeom_math::Direction2::Y),
1952            v_dom.0 - 1.0,
1953            v_dom.1 + 1.0,
1954        )?
1955        .into())
1956    };
1957    // A border this strip fitted runs along its own row or column, and its
1958    // image is that straight line. One adopted from the neighbour was fitted
1959    // at the neighbour's pace along the sweep, which is not this strip's:
1960    // its image here is read off this surface point by point.
1961    for (edge, given, straight, span) in [
1962        (&bottom, shared[0].is_some(), row_line(v_dom.0)?, u_dom),
1963        (&top, shared[1].is_some(), row_line(v_dom.1)?, u_dom),
1964        (&rail0, shared[2].is_some(), column_line(u_dom.0)?, v_dom),
1965        (&rail1, shared[3].is_some(), column_line(u_dom.1)?, v_dom),
1966    ] {
1967        let (image, range) = if given {
1968            adopted_image(model, edge, &surface_geo, tol)?
1969        } else {
1970            (straight, span)
1971        };
1972        ogeom_algo::attach_pcurve(
1973            model,
1974            edge,
1975            image,
1976            surface_id,
1977            ogeom_topo::Location::identity(),
1978            range,
1979        )?;
1980    }
1981    // The rails carry the fit's honest budget: the neighbouring strip fitted
1982    // the same transported corners independently, and the weld between them
1983    // is only as tight as both fits.
1984    for edge in [&bottom, &top, &rail0, &rail1] {
1985        model.widen(edge, ogeom_core::Tolerance::new(error)?)?;
1986    }
1987
1988    let wire = ogeom_algo::make_wire(
1989        model,
1990        &[
1991            bottom.clone(),
1992            rail1.clone(),
1993            top.reversed(),
1994            rail0.reversed(),
1995        ],
1996        tol,
1997    )?
1998    .shape;
1999    let face = ogeom_algo::make_face_on(model, surface_id, std::slice::from_ref(&wire), tol)?.shape;
2000    let mid_u = f64::midpoint(u_dom.0, u_dom.1);
2001    let mid_v = f64::midpoint(v_dom.0, v_dom.1);
2002    let s_mid = surface_geo.point_at(mid_u, mid_v, tol)?;
2003    let (du, dv) = surface_geo.d1_at(mid_u, mid_v, tol)?;
2004    let natural_out = du.cross(dv).dot(s_mid - outward_hint) >= 0.0;
2005    let face = if natural_out == !hole {
2006        face
2007    } else {
2008        face.reversed()
2009    };
2010    Ok(SkinnedStrip {
2011        face,
2012        bottom,
2013        top,
2014        rail0,
2015        rail1,
2016    })
2017}
2018
2019/// The plane every point of the rows lies in, if there is one.
2020fn plane_of_rows(rows: &[Vec<Point>], tol: Tolerances) -> Option<Plane> {
2021    let first = rows.first()?;
2022    let last = rows.last()?;
2023    let origin = *first.first()?;
2024    let across = *first.last()? - origin;
2025    let along = *last.first()? - origin;
2026    let normal = across.cross(along);
2027    if normal.magnitude() <= tol.confusion() * across.magnitude().max(along.magnitude()) {
2028        return None;
2029    }
2030    let normal = Direction::new(normal, tol).ok()?;
2031    let plane = Plane::through(origin, normal);
2032    rows.iter()
2033        .flatten()
2034        .all(|p| plane.distance_to(*p) <= tol.confusion())
2035        .then_some(plane)
2036}
2037
2038/// A strip on its own exact plane: the borders fitted through the rows and
2039/// the end columns, the face on the plane.
2040#[allow(clippy::too_many_arguments, reason = "one construction, all its data")]
2041fn planar_strip(
2042    model: &mut Model,
2043    rows: &[Vec<Point>],
2044    plane: Plane,
2045    corners: (&Shape, &Shape, &Shape, &Shape),
2046    shared: [Option<&Shape>; 4],
2047    outward_hint: Point,
2048    hole: bool,
2049    tolerance: f64,
2050    tol: Tolerances,
2051) -> OgeomResult<SkinnedStrip> {
2052    // Splines, as every swept border is: the caps and the neighbouring
2053    // strips read them so, and a spline through collinear points is the
2054    // straight segment itself.
2055    let through = |points: &[Point]| -> OgeomResult<ogeom_geom::Curve> {
2056        let fitted = ogeom_geom::fit::fit_points(points, 3, tolerance * 0.5, tol)?;
2057        if !fitted.met {
2058            ogeom_bail!(
2059                NotDone,
2060                "a planar strip's border reached {} against a target of {tolerance}",
2061                fitted.error
2062            );
2063        }
2064        Ok(ogeom_geom::Curve::BSpline(fitted.curve))
2065    };
2066    let column = |i: usize| -> Vec<Point> { rows.iter().map(|row| row[i]).collect() };
2067    let last = rows[0].len() - 1;
2068    let (c00, c10, c01, c11) = corners;
2069    let border = |model: &mut Model,
2070                  given: Option<&Shape>,
2071                  points: Vec<Point>,
2072                  from: &Shape,
2073                  to: &Shape|
2074     -> OgeomResult<Shape> {
2075        if let Some(edge) = given {
2076            return Ok(edge.clone());
2077        }
2078        let curve = through(&points)?;
2079        let domain = curve.domain();
2080        Ok(make_edge_between(model, curve, domain, from, to, tol)?.shape)
2081    };
2082    let bottom = border(model, shared[0], rows[0].clone(), c00, c10)?;
2083    let top = border(model, shared[1], rows[rows.len() - 1].clone(), c01, c11)?;
2084    let rail0 = border(model, shared[2], column(0), c00, c01)?;
2085    let rail1 = border(model, shared[3], column(last), c10, c11)?;
2086
2087    // The loop runs across, up, back and down: counter-clockwise about the
2088    // normal from the first row's run to the first column's.
2089    let across = rows[0][last] - rows[0][0];
2090    let up = rows[rows.len() - 1][0] - rows[0][0];
2091    let normal = Direction::new(across.cross(up), tol)?;
2092    let wound = Plane::through(plane.origin(), normal);
2093    let reach = rows
2094        .iter()
2095        .flatten()
2096        .map(|p| p.distance(plane.origin()))
2097        .fold(1.0_f64, f64::max)
2098        * 2.0;
2099    let surface: SurfaceGeometry =
2100        PlaneSurface::over(wound, (-reach, reach), (-reach, reach))?.into();
2101    let face = ogeom_algo::make_face_with_pcurves(
2102        model,
2103        surface.clone(),
2104        &[vec![
2105            bottom.clone(),
2106            rail1.clone(),
2107            top.reversed(),
2108            rail0.reversed(),
2109        ]],
2110        tol,
2111    )?
2112    .shape;
2113    let mid = rows[rows.len() / 2][last / 2];
2114    let natural_out = normal.vector().dot(mid - outward_hint) >= 0.0;
2115    let face = if natural_out == !hole {
2116        face
2117    } else {
2118        face.reversed()
2119    };
2120    Ok(SkinnedStrip {
2121        face,
2122        bottom,
2123        top,
2124        rail0,
2125        rail1,
2126    })
2127}
2128
2129/// A loft through circles standing coaxial on parallel planes: the solid of
2130/// revolution of the meridian through their radii, a spline through them
2131/// in the half-plane of the first circle's start. `None` where the
2132/// sections are anything else.
2133fn coaxial_circles_loft(
2134    model: &mut Model,
2135    sections: &[Shape],
2136    tol: Tolerances,
2137) -> OgeomResult<Option<Built>> {
2138    let mut circles = Vec::with_capacity(sections.len());
2139    for wire in sections {
2140        if model.kind_of(wire)? != ShapeType::Wire {
2141            return Ok(None);
2142        }
2143        let edges = model.ordered_children_of(wire)?;
2144        let [edge] = edges.as_slice() else {
2145            return Ok(None);
2146        };
2147        let (curve, _) = spine_curve_of(model, edge)?;
2148        let ogeom_geom::Curve::Circle(c) = curve else {
2149            return Ok(None);
2150        };
2151        let placed = c
2152            .circle()
2153            .transformed(&edge.transform(model.datums())?, tol)?;
2154        circles.push(placed);
2155    }
2156    let first = circles[0].frame();
2157    let (c0, z0) = (first.origin(), first.z().vector());
2158    let last = circles[circles.len() - 1].centre();
2159    let rise = last - c0;
2160    if rise.magnitude() <= tol.confusion() {
2161        return Ok(None);
2162    }
2163    let z = rise / rise.magnitude();
2164    if z.cross(z0).magnitude() > tol.angular() {
2165        return Ok(None);
2166    }
2167    let mut heights = Vec::with_capacity(circles.len());
2168    for c in &circles {
2169        let off = c.centre() - c0;
2170        if off.cross(z).magnitude() > tol.confusion() * 10.0
2171            || c.frame().z().vector().cross(z).magnitude() > tol.angular()
2172        {
2173            return Ok(None);
2174        }
2175        heights.push(off.dot(z));
2176    }
2177    if heights.windows(2).any(|w| w[1] <= w[0] + tol.confusion()) {
2178        return Ok(None);
2179    }
2180    let x = first.x().vector();
2181    let meridian: Vec<Point> = circles
2182        .iter()
2183        .zip(&heights)
2184        .map(|(c, h)| c0 + z * *h + x * c.radius())
2185        .collect();
2186    let degree = (meridian.len() - 1).min(3);
2187    let fitted = ogeom_geom::fit::fit_points(&meridian, degree, tol.confusion() * 1e-3, tol)?;
2188    let spline: ogeom_geom::Curve = fitted.curve.into();
2189    let domain = spline.domain();
2190    let top = c0 + z * heights[heights.len() - 1];
2191    let vertex = |model: &mut Model, p: Point| ogeom_algo::make_vertex(model, p).shape;
2192    let (v_axis0, v_axis1) = (vertex(model, c0), vertex(model, top));
2193    let (v_rim0, v_rim1) = (
2194        vertex(model, meridian[0]),
2195        vertex(model, meridian[meridian.len() - 1]),
2196    );
2197    let segment =
2198        |model: &mut Model, a: (&Shape, Point), b: (&Shape, Point)| -> OgeomResult<Shape> {
2199            let line: ogeom_geom::Curve = LineCurve::segment(a.1, b.1, tol)?.into();
2200            let range = line.domain();
2201            Ok(make_edge_between(model, line, range, a.0, b.0, tol)?.shape)
2202        };
2203    let bottom = segment(model, (&v_axis0, c0), (&v_rim0, meridian[0]))?;
2204    let side = make_edge_between(model, spline, domain, &v_rim0, &v_rim1, tol)?.shape;
2205    let top_edge = segment(
2206        model,
2207        (&v_rim1, meridian[meridian.len() - 1]),
2208        (&v_axis1, top),
2209    )?;
2210    let axis_edge = segment(model, (&v_axis1, top), (&v_axis0, c0))?;
2211    let wire = ogeom_algo::make_wire(model, &[bottom, side, top_edge, axis_edge], tol)?.shape;
2212    // Framed from a point inside the profile: a plane's own origin is
2213    // where its face is read when it carries no trims.
2214    let inside = c0
2215        + z * (heights[heights.len() - 1] * 0.5)
2216        + x * (circles
2217            .iter()
2218            .map(|c| c.radius())
2219            .fold(f64::INFINITY, f64::min)
2220            * 0.5);
2221    let plane = Plane::new(Frame::new(
2222        inside,
2223        Direction::new(z.cross(x), tol)?,
2224        Direction::new(x, tol)?,
2225        tol,
2226    )?);
2227    let face = ogeom_algo::make_face(model, PlaneSurface::new(plane).into(), &[wire], tol)?.shape;
2228    let axis = ogeom_math::Axis {
2229        location: c0,
2230        direction: Direction::new(z, tol)?,
2231    };
2232    let built = ogeom_algo::make_revolution(model, &face, axis, core::f64::consts::TAU, tol)?;
2233    Ok(Some(built))
2234}
2235
2236/// A loft through sections of one edge count, every vertex a corner: one
2237/// strip per edge through all the sections, meeting its neighbours along
2238/// seams through the matched corners, each strip a plane wherever its rows
2239/// share one. `None` where the sections do not pair edge for edge.
2240fn cornered_loft(
2241    model: &mut Model,
2242    sections: &[Shape],
2243    tolerance: f64,
2244    tol: Tolerances,
2245) -> OgeomResult<Option<Built>> {
2246    let mut rings: Vec<Vec<Shape>> = Vec::with_capacity(sections.len());
2247    for wire in sections {
2248        if model.kind_of(wire)? != ShapeType::Wire || !ogeom_algo::is_wire_closed(model, wire, tol)?
2249        {
2250            return Ok(None);
2251        }
2252        rings.push(model.ordered_children_of(wire)?);
2253    }
2254    let count = rings[0].len();
2255    if count < 2 || rings.iter().any(|r| r.len() != count) {
2256        return Ok(None);
2257    }
2258    let (Some(plane0), Some(plane1)) = (
2259        ogeom_algo::find_plane(model, &sections[0], tol)?,
2260        ogeom_algo::find_plane(model, &sections[sections.len() - 1], tol)?,
2261    ) else {
2262        return Ok(None);
2263    };
2264    const ALONG: usize = 16;
2265    // Every edge of every section sampled in its ring's sense.
2266    let mut samples: Vec<Vec<Vec<Point>>> = Vec::with_capacity(rings.len());
2267    for ring in &rings {
2268        let mut per_edge = Vec::with_capacity(count);
2269        for edge in ring {
2270            let (curve, range) = spine_curve_of(model, edge)?;
2271            let curve = curve.transformed(&edge.transform(model.datums())?, tol)?;
2272            let reversed = edge.orientation() == ogeom_topo::Orientation::Reversed;
2273            let mut row = Vec::with_capacity(ALONG + 1);
2274            for k in 0..=ALONG {
2275                #[allow(clippy::cast_precision_loss)]
2276                let f = k as f64 / ALONG as f64;
2277                let t = if reversed {
2278                    range.1 - (range.1 - range.0) * f
2279                } else {
2280                    range.0 + (range.1 - range.0) * f
2281                };
2282                row.push(curve.point_at(t, tol)?);
2283            }
2284            per_edge.push(row);
2285        }
2286        samples.push(per_edge);
2287    }
2288    let corners = |model: &mut Model, s: usize| -> Vec<Shape> {
2289        (0..count)
2290            .map(|e| ogeom_algo::make_vertex(model, samples[s][e][0]).shape)
2291            .collect()
2292    };
2293    let (from, to) = (corners(model, 0), corners(model, sections.len() - 1));
2294    let middle = &samples[sections.len() / 2];
2295    let hint = {
2296        let all: Vec<Point> = middle.iter().flatten().copied().collect();
2297        #[allow(clippy::cast_precision_loss)]
2298        let n = all.len() as f64;
2299        Point::from_vector(
2300            all.iter()
2301                .fold(Vector::new(0.0, 0.0, 0.0), |acc, p| acc + p.to_vector())
2302                / n,
2303        )
2304    };
2305    let mut faces = Vec::with_capacity(count + 2);
2306    let (mut bottoms, mut tops) = (Vec::with_capacity(count), Vec::with_capacity(count));
2307    let mut first_rail: Option<Shape> = None;
2308    let mut prev_rail: Option<Shape> = None;
2309    for e in 0..count {
2310        let rows: Vec<Vec<Point>> = samples.iter().map(|s| s[e].clone()).collect();
2311        let next = (e + 1) % count;
2312        let last_rail = if e + 1 == count {
2313            first_rail.clone()
2314        } else {
2315            None
2316        };
2317        let strip = skinned_strip(
2318            model,
2319            &rows,
2320            (&from[e], &from[next], &to[e], &to[next]),
2321            [None, None, prev_rail.as_ref(), last_rail.as_ref()],
2322            hint,
2323            false,
2324            tolerance,
2325            tol,
2326        )?;
2327        if e == 0 {
2328            first_rail = Some(strip.rail0.clone());
2329        }
2330        prev_rail = Some(strip.rail1.clone());
2331        faces.push(strip.face.clone());
2332        bottoms.push(strip.bottom);
2333        tops.push(strip.top);
2334    }
2335    let towards = hint - samples[0][0][0];
2336    let n0 = plane0.normal().vector();
2337    let n0 = if n0.dot(towards) > 0.0 { -n0 } else { n0 };
2338    let away = hint - samples[sections.len() - 1][0][0];
2339    let n1 = plane1.normal().vector();
2340    let n1 = if n1.dot(away) > 0.0 { -n1 } else { n1 };
2341    faces.push(plane_cap(model, samples[0][0][0], n0, &[bottoms], tol)?);
2342    faces.push(plane_cap(
2343        model,
2344        samples[sections.len() - 1][0][0],
2345        n1,
2346        &[tops],
2347        tol,
2348    )?);
2349    let sewn = sew(model, &faces, tol)?;
2350    if sewn.shells.len() != 1 || !ogeom_algo::is_shell_closed(model, &sewn.shells[0])? {
2351        ogeom_bail!(Construction, "the cornered loft did not close");
2352    }
2353    Ok(Some(make_solid(model, &sewn.shells)?))
2354}
2355
2356/// A planar cap through `at`, facing `outward`, bounded by loops of spline
2357/// or line edges lying in it; each edge's trim is its exact projection into
2358/// the plane's chart.
2359fn plane_cap(
2360    model: &mut Model,
2361    at: Point,
2362    outward: Vector,
2363    loops: &[Vec<Shape>],
2364    tol: Tolerances,
2365) -> OgeomResult<Shape> {
2366    let cap_plane = Plane::through(at, Direction::new(outward, tol)?);
2367    let mut reach = 1.0_f64;
2368    for edges in loops {
2369        for edge in edges {
2370            let (curve, range) = spine_curve_of(model, edge)?;
2371            for k in 0..8 {
2372                let p = curve.point_at(range.0 + (range.1 - range.0) * f64::from(k) / 8.0, tol)?;
2373                reach = reach.max(p.distance(at) * 2.0);
2374            }
2375        }
2376    }
2377    let surface: SurfaceGeometry =
2378        PlaneSurface::over(cap_plane, (-reach, reach), (-reach, reach))?.into();
2379    let mut wires = Vec::with_capacity(loops.len());
2380    for edges in loops {
2381        wires.push(ogeom_algo::make_wire(model, edges, tol)?.shape);
2382    }
2383    let face = ogeom_algo::make_face(model, surface, &wires, tol)?.shape;
2384    let cap_id = {
2385        let Some(ogeom_topo::NodeData::Face(data)) = model.node(&face).map(|n| n.data()) else {
2386            ogeom_bail!(Construction, "the cap holds no face data");
2387        };
2388        data.surface
2389    };
2390    let frame = cap_plane.frame();
2391    let flat = |p: Point| {
2392        let local = frame.to_local(p);
2393        Point2::new(local.x, local.y)
2394    };
2395    for edges in loops {
2396        for edge in edges {
2397            let (curve, range) = spine_curve_of(model, edge)?;
2398            let pcurve: ogeom_geom::PlanarCurve = match &curve {
2399                ogeom_geom::Curve::BSpline(bs) => {
2400                    let control2: Vec<Point2> = bs
2401                        .control_points()
2402                        .iter()
2403                        .map(|w| flat(w.point()))
2404                        .collect();
2405                    ogeom_geom::BSpline2d::new(bs.knots().clone(), control2, tol)?.into()
2406                }
2407                ogeom_geom::Curve::Line(line) => {
2408                    let axis = line.axis();
2409                    let origin = flat(axis.location);
2410                    let ahead = flat(axis.location + axis.direction.vector());
2411                    ogeom_geom::Line2d::over(
2412                        ogeom_math::Axis2::through(origin, ahead, tol)?,
2413                        range.0,
2414                        range.1,
2415                    )?
2416                    .into()
2417                }
2418                _ => ogeom_bail!(Construction, "a cap edge is neither a spline nor a line"),
2419            };
2420            ogeom_algo::attach_pcurve(
2421                model,
2422                edge,
2423                pcurve,
2424                cap_id,
2425                ogeom_topo::Location::identity(),
2426                range,
2427            )?;
2428        }
2429    }
2430    Ok(face)
2431}
2432
2433/// Loft a solid through many closed planar sections, skinned smoothly.
2434///
2435/// The sections are sampled at matched arc-length fractions from their own
2436/// traversal starts (aligning those starts is the caller's authorship),
2437/// and the skin holds every section to `tolerance`. The caps are the first
2438/// and last sections' own planes.
2439///
2440/// # Errors
2441///
2442/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if fewer than
2443/// two sections, a section is open or not planar;
2444/// [`OgeomError::NotDone`](ogeom_core::OgeomError::NotDone) if
2445/// the skin cannot reach the tolerance.
2446pub fn make_loft_skinned(
2447    model: &mut Model,
2448    sections: &[Shape],
2449    tolerance: f64,
2450    tol: Tolerances,
2451) -> OgeomResult<Built> {
2452    if sections.len() < 2 {
2453        ogeom_bail!(Construction, "a loft needs at least two sections");
2454    }
2455    let to_point = model.kind_of(&sections[sections.len() - 1])? == ShapeType::Vertex;
2456    // The smooth skin through two sections is the ruled one, and that is
2457    // built exactly: a drum or cone between circles, planes between
2458    // polygons.
2459    if sections.len() == 2
2460        && !to_point
2461        && let Ok(built) = make_loft(model, &sections[0], &sections[1], tol)
2462    {
2463        return Ok(built);
2464    }
2465    if sections.len() > 2 && !to_point {
2466        if let Some(built) = coaxial_circles_loft(model, sections, tol)? {
2467            return Ok(built);
2468        }
2469        if let Some(built) = cornered_loft(model, sections, tolerance, tol)? {
2470            return Ok(built);
2471        }
2472    }
2473    const AROUND: usize = 48;
2474    // A trailing vertex is the apex form: the skin narrows to a point and
2475    // the solid closes there without a cap.
2476    let apex = match model.kind_of(&sections[sections.len() - 1])? {
2477        ShapeType::Vertex => {
2478            if sections.len() < 2 {
2479                ogeom_bail!(
2480                    Construction,
2481                    "a loft to a point needs a section to start from"
2482                );
2483            }
2484            let Some(data) = model
2485                .node(&sections[sections.len() - 1])
2486                .and_then(|n| n.data().as_vertex())
2487            else {
2488                ogeom_bail!(Construction, "the apex vertex holds no point");
2489            };
2490            Some(data.point)
2491        }
2492        _ => None,
2493    };
2494    let wires = &sections[..sections.len() - usize::from(apex.is_some())];
2495    let mut rows: Vec<Vec<Point>> = Vec::with_capacity(sections.len());
2496    let mut cap_planes: Vec<Option<Plane>> = Vec::with_capacity(wires.len());
2497    for wire in wires {
2498        if model.kind_of(wire)? != ShapeType::Wire {
2499            ogeom_bail!(Construction, "a loft section is a wire");
2500        }
2501        if !ogeom_algo::is_wire_closed(model, wire, tol)? {
2502            ogeom_bail!(Construction, "a loft section must be closed");
2503        }
2504        // Planarity is a *cap's* requirement, not the fit's: only the
2505        // sections a cap will stand on must hold a plane. A wavy middle
2506        // section skins fine.
2507        cap_planes.push(ogeom_algo::find_plane(model, wire, tol)?);
2508        rows.push(sample_wire(model, wire, AROUND, tol)?);
2509    }
2510    // A planar end is capped by its plane; one that is not (a wavy rim),
2511    // by a patch skinned from the ring to a point inside it.
2512    let outward_at = |rows: &[Vec<Point>], planes: &[Option<Plane>], end: bool| -> EndCap {
2513        let (i, j) = if end {
2514            (rows.len() - 1, rows.len() - 2)
2515        } else {
2516            (0, 1)
2517        };
2518        let Some(plane) = &planes[i] else {
2519            return EndCap::Skinned;
2520        };
2521        let towards = rows[j][0] - rows[i][0];
2522        let n = plane.normal().vector();
2523        EndCap::Plane(if n.dot(towards) > 0.0 { -n } else { n })
2524    };
2525    let mut built = if let Some(apex) = apex {
2526        if rows.len() < 2 {
2527            // One ring to a point is exact machinery's job when it can be;
2528            // the skin still needs two rows to shape the wall, so a middle
2529            // row is interpolated halfway toward the apex.
2530            let half: Vec<Point> = rows[0]
2531                .iter()
2532                .map(|p| Point::from_vector((p.to_vector() + apex.to_vector()) * 0.5))
2533                .collect();
2534            rows.push(half);
2535        }
2536        let outward0 = match outward_at(&rows, &cap_planes, false) {
2537            EndCap::Plane(n) => n,
2538            EndCap::Skinned => {
2539                ogeom_bail!(
2540                    Construction,
2541                    "a loft to a point starts from a planar section; a cap stands on it"
2542                );
2543            }
2544        };
2545        rows.push(vec![apex; AROUND]);
2546        skinned_solid_to_apex(model, &rows, outward0, tolerance, tol)?
2547    } else {
2548        let outward0 = outward_at(&rows, &cap_planes, false);
2549        let outward1 = outward_at(&rows, &cap_planes, true);
2550        skinned_solid(model, &rows, (outward0, outward1), tolerance, tol)?
2551    };
2552    for section in sections {
2553        built.history.generate(section, built.shape.clone());
2554    }
2555    Ok(built)
2556}
2557
2558/// Sample a closed wire at `count` matched arc-length fractions.
2559fn sample_wire(
2560    model: &Model,
2561    wire: &Shape,
2562    count: usize,
2563    tol: Tolerances,
2564) -> OgeomResult<Vec<Point>> {
2565    sample_wire_from(model, wire, count, None, tol)
2566}
2567
2568/// As [`sample_wire`], with the arc-length origin rotated to the dense
2569/// sample nearest `start_hint`: how a caller says which point of each
2570/// section rows up with which, instead of leaning on traversal starts.
2571fn sample_wire_from(
2572    model: &Model,
2573    wire: &Shape,
2574    count: usize,
2575    start_hint: Option<Point>,
2576    tol: Tolerances,
2577) -> OgeomResult<Vec<Point>> {
2578    // Dense polyline by traversal, then resample by cumulative length.
2579    let mut dense: Vec<Point> = Vec::new();
2580    for edge in explore(model, wire, Filter::OfType(ShapeType::Edge))? {
2581        let Some(data) = model.node(&edge).and_then(|n| n.data().as_edge()) else {
2582            ogeom_bail!(Construction, "a section edge holds no data");
2583        };
2584        let Some(EdgeRepr::Curve3d { curve, range, .. }) = data.curve3d() else {
2585            ogeom_bail!(Construction, "a section edge has no curve");
2586        };
2587        let Some(geometry) = model.geometry().curve(*curve) else {
2588            ogeom_bail!(Dangling, "curve is not in this model");
2589        };
2590        let reversed = edge.orientation() == ogeom_topo::Orientation::Reversed;
2591        for i in 0..64 {
2592            let f = f64::from(i) / 64.0;
2593            let t = if reversed {
2594                range.1 - (range.1 - range.0) * f
2595            } else {
2596                range.0 + (range.1 - range.0) * f
2597            };
2598            dense.push(geometry.point_at(t, tol)?);
2599        }
2600    }
2601    if let Some(hint) = start_hint {
2602        let mut best = 0usize;
2603        let mut held = f64::INFINITY;
2604        for (i, p) in dense.iter().enumerate() {
2605            let d = p.distance(hint);
2606            if d < held {
2607                held = d;
2608                best = i;
2609            }
2610        }
2611        dense.rotate_left(best);
2612    }
2613    let mut lengths = vec![0.0];
2614    for w in dense.windows(2) {
2615        let last = lengths[lengths.len() - 1];
2616        lengths.push(last + w[0].distance(w[1]));
2617    }
2618    let closing = dense[dense.len() - 1].distance(dense[0]);
2619    let total = lengths[lengths.len() - 1] + closing;
2620    let mut out = Vec::with_capacity(count);
2621    let mut cursor = 0usize;
2622    for s in 0..count {
2623        #[allow(clippy::cast_precision_loss)]
2624        let target = total * (s as f64) / (count as f64);
2625        while cursor + 1 < lengths.len() && lengths[cursor + 1] < target {
2626            cursor += 1;
2627        }
2628        let (a, b) = (dense[cursor], dense[(cursor + 1) % dense.len()]);
2629        let la = lengths[cursor];
2630        let lb = if cursor + 1 < lengths.len() {
2631            lengths[cursor + 1]
2632        } else {
2633            total
2634        };
2635        let f = if lb > la {
2636            (target - la) / (lb - la)
2637        } else {
2638            0.0
2639        };
2640        out.push(a + (b - a) * f.clamp(0.0, 1.0));
2641    }
2642    Ok(out)
2643}
2644
2645/// Sweep a circular profile along a free-form spine, skinned.
2646///
2647/// Frames along the spine are rotation-minimizing (the double-reflection
2648/// construction), so the tube neither twists nor kinks where the spine
2649/// bends; the skin holds the sampled circles to `tolerance`, and the caps
2650/// sit perpendicular to the spine's ends.
2651///
2652/// # Errors
2653///
2654/// As [`make_pipe`], plus [`OgeomError::NotDone`](ogeom_core::OgeomError::NotDone) if
2655/// the skin cannot reach the
2656/// tolerance.
2657pub fn make_pipe_skinned(
2658    model: &mut Model,
2659    spine: &Shape,
2660    radius: f64,
2661    tolerance: f64,
2662    tol: Tolerances,
2663) -> OgeomResult<Built> {
2664    if !radius.is_finite() || radius <= tol.confusion() {
2665        ogeom_bail!(Construction, "a pipe of radius {radius} holds nothing");
2666    }
2667    let (curve, range) = {
2668        let Some(data) = model.node(spine).and_then(|n| n.data().as_edge()) else {
2669            ogeom_bail!(Construction, "a pipe runs along an edge");
2670        };
2671        let Some(EdgeRepr::Curve3d { curve, range, .. }) = data.curve3d() else {
2672            ogeom_bail!(Construction, "the spine has no curve");
2673        };
2674        let Some(geometry) = model.geometry().curve(*curve) else {
2675            ogeom_bail!(Dangling, "curve is not in this model");
2676        };
2677        (geometry.clone(), *range)
2678    };
2679    const STATIONS: usize = 33;
2680    const AROUND: usize = 40;
2681    let mut stations: Vec<SpineStation> = Vec::with_capacity(STATIONS);
2682    for i in 0..STATIONS {
2683        #[allow(clippy::cast_precision_loss)]
2684        let t = range.0 + (range.1 - range.0) * (i as f64) / ((STATIONS - 1) as f64);
2685        let p = curve.point_at(t, tol)?;
2686        let d = curve.d1_at(t, tol)?;
2687        let m = d.magnitude();
2688        if m <= tol.confusion() {
2689            ogeom_bail!(Construction, "the spine is degenerate at {t}");
2690        }
2691        stations.push(SpineStation {
2692            at: p,
2693            tangent: d / m,
2694            edge: 0,
2695            t,
2696        });
2697    }
2698    let normals = rmf_normals(&stations);
2699    let mut rows: Vec<Vec<Point>> = Vec::with_capacity(STATIONS);
2700    for (i, station) in stations.iter().enumerate() {
2701        let x = normals[i];
2702        let y = station.tangent.cross(x);
2703        let mut row = Vec::with_capacity(AROUND);
2704        for a in 0..AROUND {
2705            #[allow(clippy::cast_precision_loss)]
2706            let ang = core::f64::consts::TAU * (a as f64) / (AROUND as f64);
2707            row.push(station.at + (x * ang.cos() + y * ang.sin()) * radius);
2708        }
2709        rows.push(row);
2710    }
2711    let mut built = skinned_solid(
2712        model,
2713        &rows,
2714        (
2715            EndCap::Plane(-stations[0].tangent),
2716            EndCap::Plane(stations[STATIONS - 1].tangent),
2717        ),
2718        tolerance,
2719        tol,
2720    )?;
2721    built.history.generate(spine, built.shape.clone());
2722    Ok(built)
2723}
2724
2725/// One sampled spine station: where the spine is and which way it runs.
2726#[derive(Clone, Copy)]
2727struct SpineStation {
2728    at: Point,
2729    /// The unit tangent, in the direction of travel.
2730    tangent: Vector,
2731    /// The spine edge this station stands on, by position in the spine.
2732    edge: usize,
2733    /// The station's parameter on that edge's curve.
2734    t: f64,
2735}
2736
2737/// One profile wire's closed shell round the spine: smooth wires skin as a
2738/// single closed face, faceted ones as one ring strip per facet.
2739#[allow(clippy::too_many_arguments, reason = "one frame, spelled out")]
2740fn closed_loop_shell(
2741    model: &mut Model,
2742    profile_loop: &Shape,
2743    edges: &[Shape],
2744    smooth: bool,
2745    stations: &[SpineStation],
2746    normals: &[Vector],
2747    frame0: (Point, Vector),
2748    tolerance: f64,
2749    tol: Tolerances,
2750) -> OgeomResult<Shape> {
2751    const AROUND: usize = 40;
2752    let (origin, x0) = frame0;
2753    let t0 = stations[0].tangent;
2754    let y0 = t0.cross(x0);
2755    if !smooth {
2756        // A faceted profile: one ring strip per profile edge; a fit cannot
2757        // speak a corner, so each facet gets its own v-closed skin and the
2758        // strips sew along the corner loops they share within tolerance.
2759        const ALONG_EDGE: usize = 8;
2760        // Outward for a ring strip means away from the spine's own line,
2761        // not from the loop's centroid: a ring's inner side *faces* the
2762        // centroid. The hint is the station the strip's midpoint rides.
2763        let mid_station = stations[stations.len() / 2].at;
2764        let mut faces = Vec::with_capacity(edges.len());
2765        // Each corner loop is one rail edge shared by the two strips that
2766        // meet along it, the wrap included.
2767        let mut rails: Vec<Option<Shape>> = vec![None; edges.len()];
2768        for (index, edge) in edges.iter().enumerate() {
2769            let (curve, range) = spine_curve_of(model, edge)?;
2770            let reversed = edge.orientation() == ogeom_topo::Orientation::Reversed;
2771            let mut flat_row: Vec<(f64, f64)> = Vec::with_capacity(ALONG_EDGE + 1);
2772            for kk in 0..=ALONG_EDGE {
2773                #[allow(clippy::cast_precision_loss)]
2774                let f = (kk as f64) / (ALONG_EDGE as f64);
2775                let t = if reversed {
2776                    range.1 - (range.1 - range.0) * f
2777                } else {
2778                    range.0 + (range.1 - range.0) * f
2779                };
2780                let p = curve.point_at(t, tol)?;
2781                flat_row.push(((p - origin).dot(x0), (p - origin).dot(y0)));
2782            }
2783            // rows[j = station][i = across the facet].
2784            let rows: Vec<Vec<Point>> = stations
2785                .iter()
2786                .enumerate()
2787                .map(|(i, station)| {
2788                    let x = normals[i];
2789                    let y = station.tangent.cross(x);
2790                    flat_row
2791                        .iter()
2792                        .map(|(a, b)| station.at + x * *a + y * *b)
2793                        .collect()
2794                })
2795                .collect();
2796            let next = (index + 1) % edges.len();
2797            let shared = [rails[index].clone(), rails[next].clone()];
2798            let (face, rail0, rail1) = skinned_ring_strip(
2799                model,
2800                &rows,
2801                mid_station,
2802                [shared[0].as_ref(), shared[1].as_ref()],
2803                tolerance,
2804                tol,
2805            )?;
2806            rails[index] = Some(rail0);
2807            rails[next] = Some(rail1);
2808            faces.push(face);
2809        }
2810        let sewn = sew(model, &faces, tol)?;
2811        if sewn.shells.len() != 1 || !ogeom_algo::is_shell_closed(model, &sewn.shells[0])? {
2812            if std::env::var_os("OGEOM_DEBUG_RING").is_some() {
2813                use ogeom_geom::Curve3d as _;
2814                eprintln!(
2815                    "RING: {} shells from {} strips",
2816                    sewn.shells.len(),
2817                    faces.len()
2818                );
2819                for shell in &sewn.shells {
2820                    for edge in ogeom_topo::explore_unique(model, shell, ShapeType::Edge)? {
2821                        let mut uses = 0;
2822                        for f in explore(model, shell, Filter::OfType(ShapeType::Face))? {
2823                            for w in model.children_of(&f)? {
2824                                for e in model.children_of(&w)? {
2825                                    if e.node() == edge.node() {
2826                                        uses += 1;
2827                                    }
2828                                }
2829                            }
2830                        }
2831                        if uses == 1
2832                            && let Some(d) = model.node(&edge).and_then(|n| n.data().as_edge())
2833                            && let Some(ogeom_topo::EdgeRepr::Curve3d { curve, range, .. }) =
2834                                d.curve3d()
2835                            && let Some(g) = model.geometry().curve(*curve)
2836                        {
2837                            eprintln!(
2838                                "RING open edge tol {:.2e}: {:?} -> {:?}",
2839                                d.tolerance.get(),
2840                                g.point_at(range.0, tol)?,
2841                                g.point_at(range.1, tol)?
2842                            );
2843                        }
2844                    }
2845                }
2846            }
2847            ogeom_bail!(Construction, "the faceted ring did not close");
2848        }
2849        return Ok(sewn.shells[0].clone());
2850    }
2851    let samples = sample_wire(model, profile_loop, AROUND, tol)?;
2852    let flat: Vec<(f64, f64)> = samples
2853        .iter()
2854        .map(|p| ((*p - origin).dot(x0), (*p - origin).dot(y0)))
2855        .collect();
2856    let rows: Vec<Vec<Point>> = stations
2857        .iter()
2858        .enumerate()
2859        .map(|(i, station)| {
2860            let x = normals[i];
2861            let y = station.tangent.cross(x);
2862            flat.iter()
2863                .map(|(a, b)| station.at + x * *a + y * *b)
2864                .collect()
2865        })
2866        .collect();
2867    closed_skinned_shell(model, &rows, tolerance, tol)
2868}
2869
2870/// Rotation-minimizing normals along the stations, by double reflection:
2871/// reflect in each chord's plane, then in the plane bisecting the tangents.
2872/// Self-contained (it needs only the station list) and shared by every
2873/// sweep that must not twist where its spine bends.
2874fn rmf_normals(stations: &[SpineStation]) -> Vec<Vector> {
2875    let mut normals: Vec<Vector> = Vec::with_capacity(stations.len());
2876    let t0 = stations[0].tangent;
2877    let seed = if t0.cross(ogeom_math::Vector::Z).magnitude() > 0.5 {
2878        ogeom_math::Vector::Z
2879    } else {
2880        ogeom_math::Vector::X
2881    };
2882    let n0 = {
2883        let v = seed - t0 * seed.dot(t0);
2884        v / v.magnitude()
2885    };
2886    normals.push(n0);
2887    for i in 1..stations.len() {
2888        let (p0, t0) = (stations[i - 1].at, stations[i - 1].tangent);
2889        let (p1, t1) = (stations[i].at, stations[i].tangent);
2890        let n = normals[i - 1];
2891        let v1 = p1 - p0;
2892        let c1 = v1.dot(v1);
2893        if c1 <= 1e-20 {
2894            // A corner's twin station: no travel to reflect through. The
2895            // normal is reflected across the corner's mitre plane instead;
2896            // for a vector square to the incoming tangent that is exactly
2897            // the parallel transport about the corner's own axis, and the
2898            // mirror symmetry is what lands both legs' sheared sections on
2899            // one ring. A planar corner's normal lies in the mitre plane
2900            // already and carries straight across; a skew corner's does
2901            // not, and carrying it unchanged is what left the far leg's
2902            // section off the mitre.
2903            let bisector = t0 + t1;
2904            let m = bisector.magnitude();
2905            if m <= 1e-12 {
2906                normals.push(n);
2907                continue;
2908            }
2909            let b = bisector / m;
2910            normals.push(n - b * (2.0 * n.dot(b)));
2911            continue;
2912        }
2913        normals.push(rmf_step(p0, t0, n, p1, t1));
2914    }
2915    normals
2916}
2917
2918/// One rotation-minimizing step: the normal `n0` at `(p0, t0)` carried to
2919/// `(p1, t1)` by double reflection. No travel means no change.
2920fn rmf_step(p0: Point, t0: Vector, n0: Vector, p1: Point, t1: Vector) -> Vector {
2921    let v1 = p1 - p0;
2922    let c1 = v1.dot(v1);
2923    if c1 <= 1e-20 {
2924        return n0;
2925    }
2926    let nl = n0 - v1 * (2.0 / c1 * v1.dot(n0));
2927    let tl = t0 - v1 * (2.0 / c1 * v1.dot(t0));
2928    let v2 = t1 - tl;
2929    let c2 = v2.dot(v2);
2930    let next = if c2 > 1e-20 {
2931        nl - v2 * (2.0 / c2 * v2.dot(nl))
2932    } else {
2933        nl
2934    };
2935    next / next.magnitude()
2936}
2937
2938/// A leg's generators, evaluated anywhere: the spine's own curve between
2939/// stations with the rotation-minimizing normal carried one step from the
2940/// station behind, and a straight extension past either end in the end
2941/// frame: the surface a mitre trims against. Parameters are station
2942/// indices; a unit beyond an end is one station spacing.
2943struct SpineWalk<'a> {
2944    /// Each spine edge's curve, range and whether it is travelled reversed.
2945    curves: Vec<(ogeom_geom::Curve, (f64, f64), bool)>,
2946    stations: &'a [SpineStation],
2947    normals: &'a [Vector],
2948}
2949
2950/// Where two legs' generators for one profile point meet at a corner: the
2951/// point, each leg's parameter, and how far the two generators actually
2952/// miss each other (zero when the corner turns in the plane).
2953struct CornerJoin {
2954    at: Point,
2955    s1: f64,
2956    s2: f64,
2957    gap: f64,
2958}
2959
2960impl SpineWalk<'_> {
2961    /// The spine point, unit tangent and frame normal at `s` within the run
2962    /// `(rs, re)`.
2963    fn frame_at(
2964        &self,
2965        s: f64,
2966        (rs, re): (usize, usize),
2967        tol: Tolerances,
2968    ) -> OgeomResult<(Point, Vector, Vector)> {
2969        let st = self.stations;
2970        let at = |i: usize| (st[i].at, st[i].tangent, self.normals[i]);
2971        #[allow(clippy::cast_precision_loss)]
2972        let (rsf, ref_) = (rs as f64, re as f64);
2973        if s <= rsf {
2974            let (p, t, n) = at(rs);
2975            let h = st[rs].at.distance(st[(rs + 1).min(re)].at);
2976            return Ok((p + t * ((s - rsf) * h), t, n));
2977        }
2978        if s >= ref_ {
2979            let (p, t, n) = at(re);
2980            let h = st[re].at.distance(st[re.saturating_sub(1).max(rs)].at);
2981            return Ok((p + t * ((s - ref_) * h), t, n));
2982        }
2983        #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
2984        let j = s.floor() as usize;
2985        #[allow(clippy::cast_precision_loss)]
2986        let f = s - j as f64;
2987        if f <= 0.0 {
2988            return Ok(at(j));
2989        }
2990        let (curve, range, reversed) = &self.curves[st[j + 1].edge];
2991        let t_from = if st[j].edge == st[j + 1].edge {
2992            st[j].t
2993        } else if *reversed {
2994            range.1
2995        } else {
2996            range.0
2997        };
2998        let t = t_from + (st[j + 1].t - t_from) * f;
2999        let p = curve.point_at(t, tol)?;
3000        let d = curve.d1_at(t, tol)?;
3001        let m = d.magnitude();
3002        if m <= tol.confusion() {
3003            ogeom_bail!(Construction, "the spine is degenerate at {t}");
3004        }
3005        let tangent = if *reversed { -(d / m) } else { d / m };
3006        let n = rmf_step(st[j].at, st[j].tangent, self.normals[j], p, tangent);
3007        Ok((p, tangent, n))
3008    }
3009
3010    /// The generator of profile point `(a, b)` at `s` within the run.
3011    fn generator(
3012        &self,
3013        s: f64,
3014        run: (usize, usize),
3015        (a, b): (f64, f64),
3016        tol: Tolerances,
3017    ) -> OgeomResult<Point> {
3018        let (p, t, x) = self.frame_at(s, run, tol)?;
3019        let y = t.cross(x);
3020        Ok(p + x * a + y * b)
3021    }
3022
3023    /// Where the generators of one profile point on the leg `before` and
3024    /// the leg `after` a corner meet: Gauss–Newton on both parameters from
3025    /// the corner itself, minimising the distance between the two.
3026    fn join(
3027        &self,
3028        before: (usize, usize),
3029        after: (usize, usize),
3030        ab: (f64, f64),
3031        tol: Tolerances,
3032    ) -> OgeomResult<CornerJoin> {
3033        const STEP: f64 = 1e-4;
3034        #[allow(clippy::cast_precision_loss)]
3035        let (mut s1, mut s2) = (before.1 as f64, after.0 as f64);
3036        for _ in 0..60 {
3037            let g1 = self.generator(s1, before, ab, tol)?;
3038            let g2 = self.generator(s2, after, ab, tol)?;
3039            let f = g1 - g2;
3040            let d1 = (self.generator(s1 + STEP, before, ab, tol)?
3041                - self.generator(s1 - STEP, before, ab, tol)?)
3042                / (2.0 * STEP);
3043            let d2 = (self.generator(s2 + STEP, after, ab, tol)?
3044                - self.generator(s2 - STEP, after, ab, tol)?)
3045                / (2.0 * STEP);
3046            let (a11, a12, a22) = (d1.dot(d1), -d1.dot(d2), d2.dot(d2));
3047            let (b1, b2) = (-f.dot(d1), f.dot(d2));
3048            let det = a11 * a22 - a12 * a12;
3049            if det.abs() <= 1e-30 {
3050                break;
3051            }
3052            let e1 = (b1 * a22 - a12 * b2) / det;
3053            let e2 = (a11 * b2 - a12 * b1) / det;
3054            s1 += e1;
3055            s2 += e2;
3056            if e1.abs().max(e2.abs()) <= 1e-12 {
3057                break;
3058            }
3059        }
3060        let g1 = self.generator(s1, before, ab, tol)?;
3061        let g2 = self.generator(s2, after, ab, tol)?;
3062        Ok(CornerJoin {
3063            at: g1.midpoint(g2),
3064            s1,
3065            s2,
3066            gap: g1.distance(g2),
3067        })
3068    }
3069}
3070
3071/// Sweep a planar profile lying in a plane through `axis` along a helix
3072/// about `axis`: a screw motion, the profile keeping its plane through the
3073/// axis the whole way, every point of it running its own helix. The
3074/// thread and spring operation.
3075///
3076/// `pitch` is the advance per turn along `axis`, `turns` how far the
3077/// profile turns, `left_handed` turns it the other way about the axis for
3078/// the same advance, and `taper_per_turn` moves every point away from the
3079/// axis by that much per turn (a conical helix; zero for a cylindrical
3080/// one). The walls are fitted through each profile edge's exact screw
3081/// images; the caps are the profile where it starts and where it ends.
3082///
3083/// # Errors
3084///
3085/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if
3086/// the profile is not a planar face whose plane holds the axis, reaches
3087/// the axis, would meet itself one turn on (its extent along the axis is
3088/// not less than the pitch), or tapers onto the axis; if `pitch` or
3089/// `turns` is not positive.
3090/// [`OgeomError::NotDone`](ogeom_core::OgeomError::NotDone) if a wall
3091/// cannot be fitted.
3092#[allow(clippy::too_many_arguments, clippy::too_many_lines)]
3093pub fn make_helical_sweep(
3094    model: &mut Model,
3095    profile: &Shape,
3096    axis: ogeom_math::Axis,
3097    pitch: f64,
3098    turns: f64,
3099    left_handed: bool,
3100    taper_per_turn: f64,
3101    tol: Tolerances,
3102) -> OgeomResult<Built> {
3103    if !(pitch.is_finite() && pitch > 0.0) || !(turns.is_finite() && turns > 0.0) {
3104        ogeom_bail!(
3105            Construction,
3106            "a helical sweep needs a positive pitch and turn count; got {pitch} and {turns}"
3107        );
3108    }
3109    if !taper_per_turn.is_finite() {
3110        ogeom_bail!(Construction, "a taper of {taper_per_turn} is not a length");
3111    }
3112    if model.kind_of(profile)? != ShapeType::Face {
3113        ogeom_bail!(Construction, "a helical sweep sweeps a planar face");
3114    }
3115    let Some(plane) = ogeom_algo::find_plane(model, profile, tol)? else {
3116        ogeom_bail!(Construction, "a helical sweep sweeps a planar face");
3117    };
3118    let z = axis.direction.vector();
3119    if plane.normal().vector().dot(z).abs() > tol.angular()
3120        || plane.distance_to(axis.location) > tol.confusion() * 100.0
3121    {
3122        ogeom_bail!(
3123            Construction,
3124            "the profile's plane does not hold the axis; a helical sweep \
3125             turns a profile about an axis in its own plane"
3126        );
3127    }
3128    let total = core::f64::consts::TAU * turns;
3129    let sense = if left_handed { -1.0 } else { 1.0 };
3130    // The screw image of a point after turning through `theta`.
3131    let screw = |p: Point, theta: f64| -> OgeomResult<Point> {
3132        let foot = axis.project(p);
3133        let out = p - foot;
3134        let rho = out.magnitude();
3135        let grown = rho + taper_per_turn * theta / core::f64::consts::TAU;
3136        if rho <= tol.confusion() || grown <= tol.confusion() {
3137            ogeom_bail!(
3138                Construction,
3139                "the profile reaches the axis, where a helical sweep has no \
3140                 helix to follow"
3141            );
3142        }
3143        let radial = out / rho;
3144        let across = z.cross(radial);
3145        let (sin, cos) = (sense * theta).sin_cos();
3146        let turned = radial * cos + across * sin;
3147        Ok(foot + z * (pitch * theta / core::f64::consts::TAU) + turned * grown)
3148    };
3149
3150    // The profile's loops, their edges in ring order, each sampled.
3151    let loops = explore(model, profile, Filter::OfType(ShapeType::Wire))?;
3152    if loops.is_empty() {
3153        ogeom_bail!(Construction, "the profile has no loop to sweep");
3154    }
3155    // One turn on, the profile must clear itself.
3156    if turns > 1.0 {
3157        let mut low = f64::INFINITY;
3158        let mut high = f64::NEG_INFINITY;
3159        for wire in &loops {
3160            for p in sample_wire(model, wire, 64, tol)? {
3161                let h = (p - axis.location).dot(z);
3162                low = low.min(h);
3163                high = high.max(h);
3164            }
3165        }
3166        if high - low >= pitch - tol.confusion() {
3167            ogeom_bail!(
3168                Construction,
3169                "the profile spans {} along the axis, not less than the pitch \
3170                 {pitch}; a turn on it meets itself",
3171                high - low
3172            );
3173        }
3174    }
3175    let tolerance = tol.confusion() * 100.0;
3176    // The walls in quarter turns, each a strip of its own sharing its
3177    // borders with the next: one fit down many turns of a helix cannot
3178    // reach the tolerance, a quarter turn's can.
3179    #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
3180    let segments = ((turns * 4.0).ceil() as usize).max(1);
3181    // The fit keeps fewer controls than samples, so its reach at the
3182    // samples is set by how many there are.
3183    const PER_SEGMENT: usize = 48;
3184    #[allow(clippy::cast_precision_loss)]
3185    let theta_at = |seg: usize, i: usize| {
3186        total * ((seg * PER_SEGMENT + i) as f64) / ((segments * PER_SEGMENT) as f64)
3187    };
3188
3189    let mut faces: Vec<Shape> = Vec::new();
3190    let mut cap_loops: [Vec<Vec<Shape>>; 2] = [Vec::new(), Vec::new()];
3191    for (li, wire) in loops.iter().enumerate() {
3192        let hole = li != 0;
3193        let edges = model.ordered_children_of(wire)?;
3194        let centre = {
3195            let samples = sample_wire(model, wire, 32, tol)?;
3196            #[allow(clippy::cast_precision_loss)]
3197            let n = samples.len() as f64;
3198            let sum = samples
3199                .iter()
3200                .fold(Vector::new(0.0, 0.0, 0.0), |acc, p| acc + p.to_vector());
3201            Point::from_vector(sum / n)
3202        };
3203        // The ring's pieces: every edge in the ring's sense, a closed one
3204        // (a circle, the whole ring) cut in quarters so each strip's fit
3205        // spans a quarter turn round it at most.
3206        let mut pieces: Vec<(ogeom_geom::Curve, f64, f64)> = Vec::new();
3207        for edge in &edges {
3208            let (curve, range) = spine_curve_of(model, edge)?;
3209            let reversed = edge.orientation() == ogeom_topo::Orientation::Reversed;
3210            let (t0, t1) = if reversed { (range.1, range.0) } else { range };
3211            let closed =
3212                ogeom_algo::edge_vertices(model, edge)?.is_some_and(|(a, b)| a.is_same(&b));
3213            let parts = if closed { 4 } else { 1 };
3214            for k in 0..parts {
3215                #[allow(clippy::cast_precision_loss)]
3216                let (f0, f1) = (k as f64 / parts as f64, (k + 1) as f64 / parts as f64);
3217                pieces.push((curve.clone(), t0 + (t1 - t0) * f0, t0 + (t1 - t0) * f1));
3218            }
3219        }
3220        // Each piece's samples, in the ring's sense, ends on its corners
3221        // exactly.
3222        let starts: Vec<Point> = pieces
3223            .iter()
3224            .map(|(curve, a, _)| curve.point_at(*a, tol))
3225            .collect::<OgeomResult<_>>()?;
3226        let count = pieces.len();
3227        let mut rows0: Vec<Vec<Point>> = Vec::with_capacity(count);
3228        for (pi, (curve, a, b)) in pieces.iter().enumerate() {
3229            let along = if matches!(curve, ogeom_geom::Curve::Line(_)) {
3230                8
3231            } else {
3232                24
3233            };
3234            let mut row = Vec::with_capacity(along + 1);
3235            for k in 0..=along {
3236                #[allow(clippy::cast_precision_loss)]
3237                let f = (k as f64) / (along as f64);
3238                row.push(curve.point_at(a + (b - a) * f, tol)?);
3239            }
3240            row[0] = starts[pi];
3241            row[along] = starts[(pi + 1) % count];
3242            rows0.push(row);
3243        }
3244        let rows_of = |row0: &[Point], seg: usize| -> OgeomResult<Vec<Vec<Point>>> {
3245            (0..=PER_SEGMENT)
3246                .map(|i| {
3247                    row0.iter()
3248                        .map(|p| screw(*p, theta_at(seg, i)))
3249                        .collect::<OgeomResult<Vec<Point>>>()
3250                })
3251                .collect()
3252        };
3253
3254        // A vertex set at every segment boundary.
3255        let mut corners: Vec<Vec<Shape>> = Vec::with_capacity(segments + 1);
3256        for b in 0..=segments {
3257            let theta = if b == segments { total } else { theta_at(b, 0) };
3258            let mut set = Vec::with_capacity(count);
3259            for p in &starts {
3260                set.push(ogeom_algo::make_vertex(model, screw(*p, theta)?).shape);
3261            }
3262            corners.push(set);
3263        }
3264        let mut bottoms = Vec::with_capacity(count);
3265        let mut tops = Vec::with_capacity(count);
3266        // Per edge, the previous segment's top border, for the next to
3267        // start on.
3268        let mut held_tops: Vec<Option<Shape>> = vec![None; count];
3269        for seg in 0..segments {
3270            let hint = screw(centre, theta_at(seg, PER_SEGMENT / 2))?;
3271            let mut first_rail: Option<Shape> = None;
3272            let mut prev_rail: Option<Shape> = None;
3273            for ei in 0..count {
3274                let rows = rows_of(&rows0[ei], seg)?;
3275                let next = (ei + 1) % count;
3276                let last_rail = if ei + 1 == count {
3277                    first_rail.clone()
3278                } else {
3279                    None
3280                };
3281                let (from, to) = (&corners[seg], &corners[seg + 1]);
3282                let strip = skinned_strip(
3283                    model,
3284                    &rows,
3285                    (&from[ei], &from[next], &to[ei], &to[next]),
3286                    [
3287                        held_tops[ei].as_ref(),
3288                        None,
3289                        prev_rail.as_ref(),
3290                        last_rail.as_ref(),
3291                    ],
3292                    hint,
3293                    hole,
3294                    tolerance,
3295                    tol,
3296                )?;
3297                if ei == 0 {
3298                    first_rail = Some(strip.rail0.clone());
3299                }
3300                prev_rail = Some(strip.rail1.clone());
3301                faces.push(strip.face.clone());
3302                if seg == 0 {
3303                    bottoms.push(strip.bottom.clone());
3304                }
3305                if seg + 1 == segments {
3306                    tops.push(strip.top.clone());
3307                }
3308                held_tops[ei] = Some(strip.top);
3309            }
3310        }
3311        cap_loops[0].push(bottoms);
3312        cap_loops[1].push(tops);
3313    }
3314
3315    // The caps: the profile's plane where it starts, and that plane
3316    // screwed on to where it ends.
3317    for (end, loops) in cap_loops.iter().enumerate() {
3318        let theta = if end == 0 { 0.0 } else { total };
3319        let at = screw(centre_of(model, profile, tol)?, theta)?;
3320        // Square to the way the profile turns there: out of the solid,
3321        // back at the start and on at the end.
3322        let normal = {
3323            let out = at - axis.project(at);
3324            let travel = z.cross(out / out.magnitude()) * sense;
3325            if end == 0 { -travel } else { travel }
3326        };
3327        let cap_plane = Plane::through(at, Direction::new(normal, tol)?);
3328        let mut reach = 1.0_f64;
3329        for edges in loops {
3330            for edge in edges {
3331                let (curve, range) = spine_curve_of(model, edge)?;
3332                for k in 0..8 {
3333                    let p =
3334                        curve.point_at(range.0 + (range.1 - range.0) * f64::from(k) / 8.0, tol)?;
3335                    reach = reach.max(p.distance(at) * 2.0);
3336                }
3337            }
3338        }
3339        let surface: SurfaceGeometry =
3340            PlaneSurface::over(cap_plane, (-reach, reach), (-reach, reach))?.into();
3341        let mut wires = Vec::with_capacity(loops.len());
3342        for edges in loops {
3343            wires.push(ogeom_algo::make_wire(model, edges, tol)?.shape);
3344        }
3345        let face = ogeom_algo::make_face(model, surface, &wires, tol)?.shape;
3346        let cap_id = {
3347            let Some(ogeom_topo::NodeData::Face(data)) = model.node(&face).map(|n| n.data()) else {
3348                ogeom_bail!(Construction, "the cap holds no face data");
3349            };
3350            data.surface
3351        };
3352        let frame = cap_plane.frame();
3353        for edges in loops {
3354            for edge in edges {
3355                let (curve, range) = spine_curve_of(model, edge)?;
3356                let ogeom_geom::Curve::BSpline(bs) = &curve else {
3357                    ogeom_bail!(Construction, "a swept ring is not a spline");
3358                };
3359                let control2: Vec<Point2> = bs
3360                    .control_points()
3361                    .iter()
3362                    .map(|w| {
3363                        let local = frame.to_local(w.point());
3364                        Point2::new(local.x, local.y)
3365                    })
3366                    .collect();
3367                let pcurve: ogeom_geom::PlanarCurve =
3368                    ogeom_geom::BSpline2d::new(bs.knots().clone(), control2, tol)?.into();
3369                ogeom_algo::attach_pcurve(
3370                    model,
3371                    edge,
3372                    pcurve,
3373                    cap_id,
3374                    ogeom_topo::Location::identity(),
3375                    range,
3376                )?;
3377            }
3378        }
3379        faces.push(face);
3380    }
3381
3382    let sewn = sew(model, &faces, tol)?;
3383    if sewn.shells.len() != 1 || !ogeom_algo::is_shell_closed(model, &sewn.shells[0])? {
3384        ogeom_bail!(Construction, "the helical sweep did not close");
3385    }
3386    let solid = make_solid(model, &sewn.shells)?.shape;
3387    let mut history = History::new();
3388    history.generate(profile, solid.clone());
3389    Ok(Built::new(solid, history))
3390}
3391
3392/// The centroid of a face's outer ring's samples.
3393fn centre_of(model: &Model, profile: &Shape, tol: Tolerances) -> OgeomResult<Point> {
3394    let Some(wire) = explore(model, profile, Filter::OfType(ShapeType::Wire))?
3395        .into_iter()
3396        .next()
3397    else {
3398        ogeom_bail!(Construction, "the profile has no loop");
3399    };
3400    let samples = sample_wire(model, &wire, 32, tol)?;
3401    #[allow(clippy::cast_precision_loss)]
3402    let n = samples.len() as f64;
3403    let sum = samples
3404        .iter()
3405        .fold(Vector::new(0.0, 0.0, 0.0), |acc, p| acc + p.to_vector());
3406    Ok(Point::from_vector(sum / n))
3407}
3408
3409/// The pipe along a spine of lines and circular arcs meeting tangent to
3410/// one another, built exactly: down a line the section is extruded, round
3411/// an arc it is revolved about the arc's axis (the rotation-minimizing
3412/// frame of a circle is its own rotation), and the legs are fused on the
3413/// sections they share. Every wall is then the closed form its profile
3414/// edge sweeps: a plane, drum, cone, ball or torus where the edge is a
3415/// line or circle. `None` where the spine or profile is not of that kind,
3416/// for the general construction to take.
3417fn exact_legs(
3418    model: &mut Model,
3419    profile: &Shape,
3420    spine: &Shape,
3421    frenet: bool,
3422    tol: Tolerances,
3423) -> OgeomResult<Option<Built>> {
3424    use ogeom_geom::{Curve, Curve3d as _};
3425    let edges: Vec<Shape> = match model.kind_of(spine)? {
3426        ShapeType::Edge => vec![spine.clone()],
3427        ShapeType::Wire => model.ordered_children_of(spine)?,
3428        _ => return Ok(None),
3429    };
3430    let solid = match model.kind_of(profile)? {
3431        ShapeType::Face => true,
3432        ShapeType::Wire => false,
3433        _ => return Ok(None),
3434    };
3435    if edges.is_empty() || (!solid && edges.len() > 1) {
3436        return Ok(None);
3437    }
3438    // Each leg: where it starts and ends, its heading at either end, and
3439    // the motion that carries the section down it.
3440    struct Leg {
3441        start: Point,
3442        heading: (Vector, Vector),
3443        motion: Transform,
3444        along: LegKind,
3445    }
3446    enum LegKind {
3447        Line(Vector),
3448        Arc(ogeom_math::Axis, f64),
3449    }
3450    let mut legs = Vec::with_capacity(edges.len());
3451    for edge in &edges {
3452        let (curve, range) = spine_curve_of(model, edge)?;
3453        let reversed = edge.orientation() == ogeom_topo::Orientation::Reversed;
3454        let (t0, t1) = if reversed { (range.1, range.0) } else { range };
3455        let (a, b) = (curve.point_at(t0, tol)?, curve.point_at(t1, tol)?);
3456        let sense = if reversed { -1.0 } else { 1.0 };
3457        let heading = (curve.d1_at(t0, tol)? * sense, curve.d1_at(t1, tol)? * sense);
3458        let basis = match &curve {
3459            Curve::Trimmed(t) => t.basis().clone(),
3460            other => other.clone(),
3461        };
3462        let (motion, along) = match basis {
3463            Curve::Line(_) => {
3464                if frenet {
3465                    return Ok(None);
3466                }
3467                (Transform::translation(b - a), LegKind::Line(b - a))
3468            }
3469            Curve::Circle(c) => {
3470                let frame = c.circle().frame();
3471                // Turning the way the walk heads at its start.
3472                let turn = (a - frame.origin()).cross(heading.0);
3473                let direction = if turn.dot(frame.z().vector()) > 0.0 {
3474                    frame.z()
3475                } else {
3476                    -frame.z()
3477                };
3478                let axis = ogeom_math::Axis {
3479                    location: frame.origin(),
3480                    direction,
3481                };
3482                let angle = (t1 - t0).abs();
3483                (Transform::rotation(axis, angle), LegKind::Arc(axis, angle))
3484            }
3485            _ => return Ok(None),
3486        };
3487        legs.push(Leg {
3488            start: a,
3489            heading,
3490            motion,
3491            along,
3492        });
3493    }
3494    // Legs must meet tangent: a corner is mitred by the general path.
3495    for pair in legs.windows(2) {
3496        let (x, y) = (pair[0].heading.1, pair[1].heading.0);
3497        if x.cross(y).magnitude() > tol.angular() * x.magnitude() * y.magnitude() || x.dot(y) <= 0.0
3498        {
3499            return Ok(None);
3500        }
3501    }
3502    // The profile square to the spine's exact start tangent, and on it.
3503    let Some(plane) = ogeom_algo::find_plane(model, profile, tol)? else {
3504        return Ok(None);
3505    };
3506    let t0 = legs[0].heading.0;
3507    if plane.normal().vector().cross(t0).magnitude() > tol.angular() * t0.magnitude()
3508        || plane.distance_to(legs[0].start) > tol.confusion() * 100.0
3509    {
3510        return Ok(None);
3511    }
3512
3513    let mut result: Option<Shape> = None;
3514    let mut carried = Transform::IDENTITY;
3515    for leg in &legs {
3516        // Rebuilt even where it stands: the caps are this face, and it
3517        // carries its trims on its own plane.
3518        let heading = carried.apply_vector(legs[0].heading.0);
3519        let section = realized_profile_wound(model, profile, &carried, Some(heading), tol)?;
3520        let piece = match leg.along {
3521            LegKind::Line(v) => ogeom_algo::make_prism(model, &section, v, tol)?.shape,
3522            LegKind::Arc(axis, angle) => {
3523                ogeom_algo::make_revolution(model, &section, axis, angle, tol)?.shape
3524            }
3525        };
3526        carried = leg.motion * carried;
3527        result = Some(match result {
3528            None => piece,
3529            Some(held) => ogeom_bool::fuse(model, &held, &piece, tol)?.shape,
3530        });
3531    }
3532    Ok(result.map(|shape| {
3533        let mut history = History::new();
3534        history.generate(spine, shape.clone());
3535        history.generate(profile, shape.clone());
3536        Built::new(shape, history)
3537    }))
3538}
3539
3540/// Sweep a planar profile (a wire, or a face whose holes ride along)
3541/// down an arbitrary spine, one skinned wall per profile loop.
3542///
3543/// The spine may be a single edge or a wire of edges of any curve the
3544/// vocabulary evaluates: lines, arcs, splines, helices. Frames along it are
3545/// rotation-minimizing by default (the double-reflection construction), so
3546/// the profile neither twists nor kinks where the spine bends; `frenet`
3547/// asks for the Frenet frame instead, which turns with the spine's own
3548/// curvature, the law a thread wants. Stations are placed by each edge's
3549/// own turning, the skin holds every transported section to `tolerance`,
3550/// and the caps sit perpendicular to the spine's ends, holes and all.
3551///
3552/// Each spine edge skins its own run of wall, and neighbouring runs share
3553/// the section where their edges meet, so a join where the curvature steps
3554/// (an arc running on into its tangent line) is followed exactly rather
3555/// than smoothed by one fit across it. A run whose sections all lie in one
3556/// plane is that plane.
3557///
3558/// A sharp corner is mitred. Between straight legs the mitre is a plane and
3559/// each wall is sheared onto it; where a leg is curved the two legs' walls
3560/// end on the crossing of their generators (each profile point's own path
3561/// down either leg, run straight on past the corner), which is exact where
3562/// the corner turns in the leg's plane. Where it turns a curved leg out of
3563/// its plane the generators miss, and the spine is swept in pieces instead:
3564/// each side runs on straight past the corner, is trimmed by the mitre
3565/// plane, and the pieces are fused, the difference between their sections
3566/// standing as a face of the mitre plane.
3567///
3568/// # Errors
3569///
3570/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the
3571/// profile is not planar, leans along the spine, or does not sit at the
3572/// spine's start; if `frenet` is asked of a spine that never bends, or of a
3573/// cornered one; if a wire (not a face) is swept round a corner that turns
3574/// a curved leg out of its plane, which only solid pieces can mitre; if the
3575/// spine all but doubles back at a corner; or if a leg is shorter than its
3576/// corner's reach.
3577/// [`OgeomError::NotDone`](ogeom_core::OgeomError::NotDone) if the skin
3578/// cannot reach the tolerance.
3579pub fn make_pipe_shell(
3580    model: &mut Model,
3581    profile: &Shape,
3582    spine: &Shape,
3583    frenet: bool,
3584    tolerance: f64,
3585    tol: Tolerances,
3586) -> OgeomResult<Built> {
3587    const AROUND: usize = 40;
3588
3589    if let Some(exact) = exact_legs(model, profile, spine, frenet, tol)? {
3590        return Ok(exact);
3591    }
3592    let stations = shell_stations(model, spine, tol)?;
3593    // Corners: twin stations standing on one point with different headings.
3594    let corners: Vec<usize> = (0..stations.len() - 1)
3595        .filter(|&i| {
3596            stations[i].at.distance(stations[i + 1].at) <= tol.confusion()
3597                && (stations[i]
3598                    .tangent
3599                    .cross(stations[i + 1].tangent)
3600                    .magnitude()
3601                    > tol.angular()
3602                    || stations[i].tangent.dot(stations[i + 1].tangent) < 0.0)
3603        })
3604        .collect();
3605    // Every twin, corners and smooth junctions alike: each ends one run of
3606    // skin and starts the next, the two runs sharing the section there. A
3607    // smooth junction's mitre plane is its own section, so nothing shears.
3608    let kinks: Vec<usize> = (0..stations.len() - 1)
3609        .filter(|&i| stations[i].at.distance(stations[i + 1].at) <= tol.confusion())
3610        .collect();
3611    let ring = stations[0].at.distance(stations[stations.len() - 1].at) <= tol.confusion() * 10.0;
3612    if ring && kinks.is_empty() {
3613        return closed_pipe_shell(model, profile, spine, stations, frenet, tolerance, tol);
3614    }
3615    // The cornered ring closes at its own wrap. Seamed on a corner, the
3616    // wrap is one more mitre; seamed mid-leg, the wrap's "mitre" plane is
3617    // the leg's own cross-section (both twin tangents are the leg's), and
3618    // the shear onto it moves nothing, so the two halves of that leg butt
3619    // together on the seam's own ring, coplanar walls meeting on it. The
3620    // solid is exact either way; the mid-leg seam merely leaves its leg in
3621    // two pieces.
3622    if frenet && !corners.is_empty() {
3623        ogeom_bail!(
3624            Construction,
3625            "a Frenet frame has no direction at a corner; sweep a cornered \
3626             spine with the rotation-minimizing frame"
3627        );
3628    }
3629    let normals = if frenet {
3630        frenet_normals(&stations, tol)?
3631    } else {
3632        rmf_normals(&stations)
3633    };
3634    // A ring's frame must come home: carry once more across the wrap
3635    // corner, read the twist between departure and return, and spread it
3636    // along the arc: the smooth loop's own reconciliation, ending at a
3637    // mitre instead of a tangent join.
3638    let normals = if ring {
3639        let mut extended = stations.clone();
3640        extended.push(stations[0]);
3641        let carried = rmf_normals(&extended);
3642        let (n0, n_home) = (carried[0], carried[carried.len() - 1]);
3643        let t0 = stations[0].tangent;
3644        let twist = (n0.cross(n_home).dot(t0)).atan2(n0.dot(n_home));
3645        let mut lengths = vec![0.0_f64];
3646        for pair in extended.windows(2) {
3647            let held = lengths[lengths.len() - 1];
3648            lengths.push(held + pair[0].at.distance(pair[1].at));
3649        }
3650        let total = lengths[lengths.len() - 1];
3651        carried
3652            .iter()
3653            .take(stations.len())
3654            .enumerate()
3655            .map(|(i, n)| {
3656                let phi = -twist * lengths[i] / total;
3657                let t = extended[i].tangent;
3658                let v = *n * phi.cos() + t.cross(*n) * phi.sin();
3659                let v = v - t * v.dot(t);
3660                v / v.magnitude()
3661            })
3662            .collect()
3663    } else {
3664        normals
3665    };
3666    // At each corner both twin sections are thrown onto the bisector plane
3667    // along their own tangents; the mirror symmetry of the rotation-
3668    // minimizing frame lands them on one ring, the mitre both runs share.
3669    let mitre: Vec<Option<(Point, Vector)>> = {
3670        let mut out: Vec<Option<(Point, Vector)>> = vec![None; stations.len()];
3671        for &k in &kinks {
3672            let n = stations[k].tangent + stations[k + 1].tangent;
3673            if n.magnitude() <= tol.angular() {
3674                ogeom_bail!(
3675                    Construction,
3676                    "the spine doubles straight back on itself; no mitre \
3677                     plane divides that corner"
3678                );
3679            }
3680            out[k] = Some((stations[k].at, n));
3681            out[k + 1] = Some((stations[k + 1].at, n));
3682        }
3683        if ring {
3684            let wrap = stations.len() - 1;
3685            let n = stations[wrap].tangent + stations[0].tangent;
3686            if n.magnitude() <= tol.angular() {
3687                ogeom_bail!(
3688                    Construction,
3689                    "the spine doubles straight back on itself; no mitre \
3690                     plane divides that corner"
3691                );
3692            }
3693            out[wrap] = Some((stations[wrap].at, n));
3694            out[0] = Some((stations[0].at, n));
3695        }
3696        out
3697    };
3698    let runs: Vec<(usize, usize)> = {
3699        let mut out = Vec::with_capacity(kinks.len() + 1);
3700        let mut start = 0;
3701        for &k in &kinks {
3702            out.push((start, k));
3703            start = k + 1;
3704        }
3705        out.push((start, stations.len() - 1));
3706        out
3707    };
3708    // A mitred end between straight legs is a *shear*: the honest wall is
3709    // the run's own surface trimmed by the mitre plane, which for a
3710    // straight leg is exactly the ruled skin between its two end rings.
3711    let straight = |rs: usize, re: usize| -> bool {
3712        let t0 = stations[rs].tangent;
3713        (rs..=re).all(|i| stations[i].tangent.cross(t0).magnitude() <= tol.angular())
3714    };
3715    // Every corner as the pair of runs it stands between, the wrap
3716    // included, and whether either leg is curved. A curved leg's trim is
3717    // not a loft of its rows: the two legs' generators for one profile
3718    // point are followed (the leg's own curve, run straight on past the
3719    // corner) to where they meet, and each wall ends on that crossing.
3720    // Where the corner turns in the plane the crossing is exact; a skew
3721    // corner's generators miss each other, and that miss is refused.
3722    struct CornerPair {
3723        before: (usize, usize),
3724        after: (usize, usize),
3725        curved: bool,
3726    }
3727    let corner_pairs: Vec<CornerPair> = {
3728        let mut out = Vec::new();
3729        // A smooth junction joins its runs on their shared section; only a
3730        // corner that turns asks the generators where the walls meet.
3731        for pair in runs.windows(2) {
3732            out.push(CornerPair {
3733                before: pair[0],
3734                after: pair[1],
3735                curved: corners.contains(&pair[0].1)
3736                    && (!straight(pair[0].0, pair[0].1) || !straight(pair[1].0, pair[1].1)),
3737            });
3738        }
3739        if ring && runs.len() > 1 {
3740            let (before, after) = (runs[runs.len() - 1], runs[0]);
3741            let turns = stations[before.1]
3742                .tangent
3743                .cross(stations[after.0].tangent)
3744                .magnitude()
3745                > tol.angular()
3746                || stations[before.1].tangent.dot(stations[after.0].tangent) < 0.0;
3747            out.push(CornerPair {
3748                before,
3749                after,
3750                curved: turns && (!straight(before.0, before.1) || !straight(after.0, after.1)),
3751            });
3752        }
3753        out
3754    };
3755    let curves: Vec<(ogeom_geom::Curve, (f64, f64), bool)> = {
3756        let edges: Vec<Shape> = match model.kind_of(spine)? {
3757            ShapeType::Edge => vec![spine.clone()],
3758            _ => model.ordered_children_of(spine)?,
3759        };
3760        let mut out = Vec::with_capacity(edges.len());
3761        for edge in &edges {
3762            let (curve, range) = spine_curve_of(model, edge)?;
3763            out.push((
3764                curve,
3765                range,
3766                edge.orientation() == ogeom_topo::Orientation::Reversed,
3767            ));
3768        }
3769        out
3770    };
3771    let walk = SpineWalk {
3772        curves,
3773        stations: &stations,
3774        normals: &normals,
3775    };
3776    let join_reach = tolerance.max(tol.confusion() * 100.0);
3777    let curved_join = |pair: &CornerPair, ab: (f64, f64)| -> OgeomResult<CornerJoin> {
3778        let join = walk.join(pair.before, pair.after, ab, tol)?;
3779        if join.gap > join_reach {
3780            ogeom_bail!(
3781                Construction,
3782                "a skew corner against a curved leg is still owed its frame \
3783                 law: the legs' generators miss by {}; see docs/PARITY.md, \
3784                 offset.sweeps",
3785                join.gap
3786            );
3787        }
3788        Ok(join)
3789    };
3790    // Across a curved corner the two legs' skins share the join row as one
3791    // edge: a run adopts the previous run's end row at its start, and the
3792    // last run of a cornered ring adopts the first run's start row at its
3793    // end.
3794    let shares_start = |ri: usize| -> bool {
3795        ri > 0
3796            && corner_pairs
3797                .iter()
3798                .any(|pair| pair.curved && pair.after == runs[ri])
3799    };
3800    let shares_end = |ri: usize| -> bool {
3801        ring && ri + 1 == runs.len()
3802            && corner_pairs
3803                .iter()
3804                .any(|pair| pair.curved && pair.before == runs[ri] && pair.after == runs[0])
3805    };
3806    // The curved corner a station is a twin of, if any.
3807    let curved_at = |i: usize| -> Option<&CornerPair> {
3808        corner_pairs
3809            .iter()
3810            .find(|pair| pair.curved && (pair.before.1 == i || pair.after.0 == i))
3811    };
3812
3813    // The profile's loops: a face contributes every wire, holes included;
3814    // a bare wire is one loop.
3815    let loops: Vec<Shape> = match model.kind_of(profile)? {
3816        ShapeType::Face => explore(model, profile, Filter::OfType(ShapeType::Wire))?,
3817        ShapeType::Wire => vec![profile.clone()],
3818        other => ogeom_bail!(
3819            Construction,
3820            "a pipe shell sweeps a planar wire or face, not a {other:?}"
3821        ),
3822    };
3823    if loops.is_empty() {
3824        ogeom_bail!(Construction, "the profile has no loop to sweep");
3825    }
3826    let Some(plane) = ogeom_algo::find_plane(model, profile, tol)? else {
3827        ogeom_bail!(Construction, "a pipe shell sweeps a planar profile");
3828    };
3829    let t0 = stations[0].tangent;
3830    // Square to the spine's own start tangent, read exactly off its first
3831    // edge, to an angle's tolerance.
3832    let exact_t0 = {
3833        let first = match model.kind_of(spine)? {
3834            ShapeType::Edge => spine.clone(),
3835            _ => model.ordered_children_of(spine)?[0].clone(),
3836        };
3837        let (curve, range) = spine_curve_of(model, &first)?;
3838        let reversed = first.orientation() == ogeom_topo::Orientation::Reversed;
3839        let d = curve.d1_at(if reversed { range.1 } else { range.0 }, tol)?;
3840        let d = if reversed { -d } else { d };
3841        d / d.magnitude()
3842    };
3843    if plane.normal().vector().cross(exact_t0).magnitude() > tol.angular() {
3844        ogeom_bail!(
3845            Construction,
3846            "the profile leans along its spine; a pipe shell runs square to \
3847             the start"
3848        );
3849    }
3850    if plane.distance_to(stations[0].at) > tol.confusion() * 100.0 {
3851        ogeom_bail!(
3852            Construction,
3853            "the profile does not sit at the spine's start"
3854        );
3855    }
3856
3857    // Transport: each loop expressed in the start frame's own 2D
3858    // coordinates, then re-expressed in every station's frame. A loop of one
3859    // smooth closed edge skins as one wall; a faceted loop skins one strip
3860    // per edge, cornered at shared vertices, because no single fit can
3861    // speak a corner.
3862    let x0 = normals[0];
3863    let y0 = t0.cross(x0);
3864    let origin = stations[0].at;
3865    let flat = |p: Point| -> (f64, f64) { ((p - origin).dot(x0), (p - origin).dot(y0)) };
3866    // A skew corner against a curved leg: the legs' generators miss, so
3867    // no join row closes the walls. Such a corner is a mitre instead: the
3868    // spine is split there, each side swept on straight past the corner
3869    // and trimmed by the mitre plane, and the pieces fused.
3870    let skew: Vec<(usize, usize)> = {
3871        let mut probes: Vec<(f64, f64)> = Vec::new();
3872        for wire in &loops {
3873            probes.extend(sample_wire(model, wire, AROUND, tol)?.into_iter().map(flat));
3874        }
3875        let mut out = Vec::new();
3876        for pair in corner_pairs.iter().filter(|pair| pair.curved) {
3877            let mut worst = 0.0_f64;
3878            for ab in &probes {
3879                worst = worst.max(walk.join(pair.before, pair.after, *ab, tol)?.gap);
3880            }
3881            if worst > join_reach {
3882                out.push((pair.before.1, pair.after.0));
3883            }
3884        }
3885        out
3886    };
3887    if !skew.is_empty() {
3888        let probes: Vec<(f64, f64)> = {
3889            let mut out = Vec::new();
3890            for wire in &loops {
3891                out.extend(sample_wire(model, wire, AROUND, tol)?.into_iter().map(flat));
3892            }
3893            out
3894        };
3895        return mitred_pieces(
3896            model, profile, spine, &stations, &skew, ring, &probes, tolerance, tol,
3897        );
3898    }
3899    let place = |i: usize, (a, b): (f64, f64)| -> OgeomResult<Point> {
3900        if let Some(pair) = curved_at(i) {
3901            return Ok(curved_join(pair, (a, b))?.at);
3902        }
3903        let x = normals[i];
3904        let y = stations[i].tangent.cross(x);
3905        let p = stations[i].at + x * a + y * b;
3906        Ok(match mitre[i] {
3907            Some((corner, n)) => {
3908                let t = stations[i].tangent;
3909                p + t * ((corner - p).dot(n) / t.dot(n))
3910            }
3911            None => p,
3912        })
3913    };
3914    // The rows a run's skin interpolates. A straight run is its two end
3915    // rings, ruled: the trimmed prism itself, whether an end is sheared
3916    // onto a mitre plane or stands on a curved corner's crossing. A curved
3917    // run is its stations, except that a stretch at a curved corner is
3918    // re-rowed: each row runs along the generators from the last plain
3919    // station to the crossing, so the skin is the leg's own surface up to
3920    // the join and nothing past it.
3921    let run_rows =
3922        |(rs, re): (usize, usize), flat_row: &[(f64, f64)]| -> OgeomResult<Vec<Vec<Point>>> {
3923            let run = (rs, re);
3924            let joins_at = |at_start: bool| -> OgeomResult<Option<Vec<CornerJoin>>> {
3925                let pair = corner_pairs.iter().find(|pair| {
3926                    pair.curved
3927                        && if at_start {
3928                            pair.after == run
3929                        } else {
3930                            pair.before == run
3931                        }
3932                });
3933                let Some(pair) = pair else {
3934                    return Ok(None);
3935                };
3936                let mut out = Vec::with_capacity(flat_row.len());
3937                for ab in flat_row {
3938                    out.push(curved_join(pair, *ab)?);
3939                }
3940                Ok(Some(out))
3941            };
3942            let (start, end) = (joins_at(true)?, joins_at(false)?);
3943            #[allow(clippy::cast_precision_loss)]
3944            let (rsf, ref_) = (rs as f64, re as f64);
3945            let short = || {
3946                ogeom_bail!(
3947                    Construction,
3948                    "a leg is shorter than its corner's reach; the mitre would \
3949                 run off its far end"
3950                )
3951            };
3952            let mut rows: Vec<Vec<Point>> = Vec::new();
3953            if straight(rs, re) {
3954                if start
3955                    .as_ref()
3956                    .is_some_and(|js| js.iter().any(|j| j.s2 >= ref_ - 0.5))
3957                    || end
3958                        .as_ref()
3959                        .is_some_and(|js| js.iter().any(|j| j.s1 <= rsf + 0.5))
3960                {
3961                    return short();
3962                }
3963                for (at_start, joins) in [(true, &start), (false, &end)] {
3964                    rows.push(match joins {
3965                        Some(js) => js.iter().map(|j| j.at).collect(),
3966                        None => {
3967                            let i = if at_start { rs } else { re };
3968                            flat_row
3969                                .iter()
3970                                .map(|ab| place(i, *ab))
3971                                .collect::<OgeomResult<Vec<Point>>>()?
3972                        }
3973                    });
3974                }
3975                return Ok(rows);
3976            }
3977            if start.is_none() && end.is_none() {
3978                for i in rs..=re {
3979                    rows.push(
3980                        flat_row
3981                            .iter()
3982                            .map(|ab| place(i, *ab))
3983                            .collect::<OgeomResult<Vec<Point>>>()?,
3984                    );
3985                }
3986                return Ok(rows);
3987            }
3988            // A curved run with a crossing at either end is re-rowed whole:
3989            // every column runs its own generator from its start to its end,
3990            // sampled at the same fractions, so the grid's shared parameter is
3991            // honest for every column; a stretch skewed only near the corner
3992            // would pace each column differently and the fit would fight it.
3993            if start
3994                .as_ref()
3995                .is_some_and(|js| js.iter().any(|j| j.s2 >= ref_ - 0.5))
3996                || end
3997                    .as_ref()
3998                    .is_some_and(|js| js.iter().any(|j| j.s1 <= rsf + 0.5))
3999            {
4000                return short();
4001            }
4002            // Each column at equal fractions of its own arc length: the grid's
4003            // parameter is one for all columns, and a column's pace differs
4004            // between the leg's curve and its straight extension, so the
4005            // fractions are taken along the generator, not along its parameter.
4006            let steps = re - rs;
4007            let fine = steps * 8;
4008            let mut columns: Vec<Vec<Point>> = Vec::with_capacity(flat_row.len());
4009            for (k, ab) in flat_row.iter().enumerate() {
4010                let lo = start.as_ref().map_or(rsf, |js| js[k].s2);
4011                let hi = end.as_ref().map_or(ref_, |js| js[k].s1);
4012                let mut along: Vec<(f64, f64)> = Vec::with_capacity(fine + 1);
4013                let mut prev: Option<Point> = None;
4014                let mut length = 0.0;
4015                for i in 0..=fine {
4016                    #[allow(clippy::cast_precision_loss)]
4017                    let sp = lo + (hi - lo) * (i as f64) / (fine as f64);
4018                    let p = walk.generator(sp, run, *ab, tol)?;
4019                    if let Some(q) = prev {
4020                        length += q.distance(p);
4021                    }
4022                    along.push((length, sp));
4023                    prev = Some(p);
4024                }
4025                let mut column = Vec::with_capacity(steps + 1);
4026                for i in 0..=steps {
4027                    #[allow(clippy::cast_precision_loss)]
4028                    let target = length * (i as f64) / (steps as f64);
4029                    let at = along.partition_point(|(l, _)| *l < target).clamp(1, fine);
4030                    let ((l0, s0), (l1, s1)) = (along[at - 1], along[at]);
4031                    let f = if l1 > l0 {
4032                        ((target - l0) / (l1 - l0)).clamp(0.0, 1.0)
4033                    } else {
4034                        0.0
4035                    };
4036                    column.push(walk.generator(s0 + (s1 - s0) * f, run, *ab, tol)?);
4037                }
4038                columns.push(column);
4039            }
4040            for i in 0..=steps {
4041                rows.push(columns.iter().map(|c| c[i]).collect());
4042            }
4043            Ok(rows)
4044        };
4045    let last = stations.len() - 1;
4046
4047    enum LoopWall {
4048        Ring {
4049            ring0: Shape,
4050            ring1: Shape,
4051        },
4052        Chain {
4053            bottoms: Vec<Shape>,
4054            tops: Vec<Shape>,
4055        },
4056    }
4057    let mut faces: Vec<Shape> = Vec::new();
4058    let mut ends: Vec<LoopWall> = Vec::with_capacity(loops.len());
4059    for (li, wire) in loops.iter().enumerate() {
4060        let hole = li != 0;
4061        let edges = model.ordered_children_of(wire)?;
4062        let single_smooth = edges.len() == 1 && {
4063            let (curve, _) = spine_curve_of(model, &edges[0])?;
4064            !matches!(curve, ogeom_geom::Curve::Line(_))
4065                && ogeom_algo::edge_vertices(model, &edges[0])?.is_some_and(|(a, b)| a.is_same(&b))
4066        };
4067        if single_smooth {
4068            let samples = sample_wire(model, wire, AROUND, tol)?;
4069            let flat_row: Vec<(f64, f64)> = samples.iter().map(|p| flat(*p)).collect();
4070            // One wall per smooth run: a fit across a corner speaks nothing,
4071            // and the twin stations put both runs' boundary rows on the one
4072            // mitred ring, where the sew joins them.
4073            let mut ring0: Option<Shape> = None;
4074            let mut ring1: Option<Shape> = None;
4075            for (ri, &(rs, re)) in runs.iter().enumerate() {
4076                let rows = run_rows((rs, re), &flat_row)?;
4077                let shared_start = shares_start(ri).then_some(()).and(ring1.as_ref());
4078                let shared_end = shares_end(ri).then_some(()).and(ring0.as_ref());
4079                let wall = skinned_wall(model, &rows, (shared_start, shared_end), tolerance, tol)?;
4080                faces.push(if hole {
4081                    wall.face.reversed()
4082                } else {
4083                    wall.face.clone()
4084                });
4085                if ring0.is_none() {
4086                    ring0 = Some(wall.ring0);
4087                }
4088                ring1 = Some(wall.ring1);
4089            }
4090            let (Some(ring0), Some(ring1)) = (ring0, ring1) else {
4091                ogeom_bail!(Construction, "the sweep produced no wall");
4092            };
4093            ends.push(LoopWall::Ring { ring0, ring1 });
4094        } else {
4095            // Shared corner vertices at both ends of every edge junction.
4096            let count = edges.len();
4097            let mut corner_flat: Vec<(f64, f64)> = Vec::with_capacity(count);
4098            for edge in &edges {
4099                // Already in the ring's own sense: a reversed edge's
4100                // vertices come back end first.
4101                let Some((start, _)) = ogeom_algo::edge_vertices(model, edge)? else {
4102                    ogeom_bail!(Construction, "a profile edge has no vertices");
4103                };
4104                let Some(data) = model.node(&start).and_then(|n| n.data().as_vertex()) else {
4105                    ogeom_bail!(Construction, "a profile vertex holds no data");
4106                };
4107                corner_flat.push(flat(data.point));
4108            }
4109            let make_corners = |model: &mut Model, station: usize| -> OgeomResult<Vec<Shape>> {
4110                let mut out = Vec::with_capacity(corner_flat.len());
4111                for ab in &corner_flat {
4112                    out.push(ogeom_algo::make_vertex(model, place(station, *ab)?).shape);
4113                }
4114                Ok(out)
4115            };
4116            // Corner vertex sets at every run boundary; a kink's twin
4117            // stations land on the same mitred points, so both runs take
4118            // the same vertex objects.
4119            let mut corners_at: Vec<Option<Vec<Shape>>> = vec![None; stations.len()];
4120            if ring {
4121                // The wrap is one corner: both runs take the same vertex
4122                // objects. Every corner's two sheared sections must land on
4123                // one ring for the loop to close; a planar ring's do
4124                // exactly, and a skew ring's (whose parallel-carried frame
4125                // leaves the far tangent's plane) do not, so the residue
4126                // is measured and the skew ring refused by name rather
4127                // than sewn hoping.
4128                let mut worst = 0.0_f64;
4129                for ab in &corner_flat {
4130                    worst = worst.max(place(last, *ab)?.distance(place(0, *ab)?));
4131                    for &k in &kinks {
4132                        worst = worst.max(place(k, *ab)?.distance(place(k + 1, *ab)?));
4133                    }
4134                }
4135                if worst > tolerance.max(tol.confusion() * 100.0) {
4136                    ogeom_bail!(
4137                        Construction,
4138                        "a skew-cornered ring's sections do not meet on \
4139                         their mitres; the out-of-plane corner's frame law \
4140                         is still owed; see docs/PARITY.md, offset.sweeps"
4141                    );
4142                }
4143                let set = make_corners(model, 0)?;
4144                if worst > tol.confusion() {
4145                    for v in &set {
4146                        model.widen(v, ogeom_core::Tolerance::new(worst * 2.0)?)?;
4147                    }
4148                }
4149                corners_at[0] = Some(set.clone());
4150                corners_at[last] = Some(set);
4151            } else {
4152                corners_at[0] = Some(make_corners(model, 0)?);
4153                corners_at[last] = Some(make_corners(model, last)?);
4154            }
4155            for &k in &kinks {
4156                let set = make_corners(model, k)?;
4157                corners_at[k] = Some(set.clone());
4158                corners_at[k + 1] = Some(set);
4159            }
4160
4161            // The loop's own centroid line, for orienting each strip.
4162            let hint_flat = {
4163                let mut a = 0.0;
4164                let mut b = 0.0;
4165                for (fa, fb) in &corner_flat {
4166                    a += fa;
4167                    b += fb;
4168                }
4169                #[allow(clippy::cast_precision_loss)]
4170                let n = count as f64;
4171                (a / n, b / n)
4172            };
4173
4174            let mut bottoms = Vec::with_capacity(count);
4175            let mut tops = Vec::with_capacity(count);
4176            // Per run: the first strip's start rail, for the last strip to
4177            // close the loop on, and the previous strip's end rail, for
4178            // the next to start from: one edge for both, never two fits.
4179            let mut run_rails: Vec<(Option<Shape>, Option<Shape>)> = vec![(None, None); runs.len()];
4180            for (ei, edge) in edges.iter().enumerate() {
4181                let (curve, range) = spine_curve_of(model, edge)?;
4182                let reversed = edge.orientation() == ogeom_topo::Orientation::Reversed;
4183                const ALONG_EDGE: usize = 8;
4184                let mut flat_row: Vec<(f64, f64)> = Vec::with_capacity(ALONG_EDGE + 1);
4185                for k in 0..=ALONG_EDGE {
4186                    #[allow(clippy::cast_precision_loss)]
4187                    let f = (k as f64) / (ALONG_EDGE as f64);
4188                    let t = if reversed {
4189                        range.1 - (range.1 - range.0) * f
4190                    } else {
4191                        range.0 + (range.1 - range.0) * f
4192                    };
4193                    flat_row.push(flat(curve.point_at(t, tol)?));
4194                }
4195                let next = (ei + 1) % count;
4196                let mut bottom: Option<Shape> = None;
4197                let mut top: Option<Shape> = None;
4198                for (ri, &(rs, re)) in runs.iter().enumerate() {
4199                    let rows = run_rows((rs, re), &flat_row)?;
4200                    let shared_start = shares_start(ri).then_some(()).and(top.as_ref());
4201                    let shared_end = shares_end(ri).then_some(()).and(bottom.as_ref());
4202                    let mid_i = usize::midpoint(rs, re);
4203                    let hint = {
4204                        let x = normals[mid_i];
4205                        let y = stations[mid_i].tangent.cross(x);
4206                        stations[mid_i].at + x * hint_flat.0 + y * hint_flat.1
4207                    };
4208                    let (Some(from), Some(to)) = (&corners_at[rs], &corners_at[re]) else {
4209                        ogeom_bail!(Construction, "a run boundary has no corners");
4210                    };
4211                    let (first_rail0, prev_rail1) = run_rails[ri].clone();
4212                    let shared_rail0 = if ei > 0 { prev_rail1 } else { None };
4213                    let shared_rail1 = if ei + 1 == count && count > 1 {
4214                        first_rail0.clone()
4215                    } else {
4216                        None
4217                    };
4218                    let strip = skinned_strip(
4219                        model,
4220                        &rows,
4221                        (&from[ei], &from[next], &to[ei], &to[next]),
4222                        [
4223                            shared_start,
4224                            shared_end,
4225                            shared_rail0.as_ref(),
4226                            shared_rail1.as_ref(),
4227                        ],
4228                        hint,
4229                        hole,
4230                        tolerance,
4231                        tol,
4232                    )?;
4233                    run_rails[ri] = (
4234                        if ei == 0 {
4235                            Some(strip.rail0.clone())
4236                        } else {
4237                            first_rail0
4238                        },
4239                        Some(strip.rail1.clone()),
4240                    );
4241                    faces.push(strip.face.clone());
4242                    if bottom.is_none() {
4243                        bottom = Some(strip.bottom);
4244                    }
4245                    top = Some(strip.top);
4246                }
4247                let (Some(bottom), Some(top)) = (bottom, top) else {
4248                    ogeom_bail!(Construction, "the sweep produced no strip");
4249                };
4250                bottoms.push(bottom);
4251                tops.push(top);
4252            }
4253            ends.push(LoopWall::Chain { bottoms, tops });
4254        }
4255    }
4256
4257    // A cap per end: one plane, one wire per loop, each edge's pcurve the
4258    // exact projection of its control net into the plane's chart. A ring
4259    // has no ends: its two boundary rings stand on one mitre plane and the
4260    // sew joins them.
4261    for end in 0..if ring { 0 } else { 2 } {
4262        let (at, outward) = if end == 0 {
4263            (stations[0].at, -stations[0].tangent)
4264        } else {
4265            (stations[last].at, stations[last].tangent)
4266        };
4267        let cap_plane = Plane::through(at, Direction::new(outward, tol)?);
4268        let mut loop_edges: Vec<Vec<Shape>> = Vec::with_capacity(ends.len());
4269        for wall in &ends {
4270            loop_edges.push(match wall {
4271                LoopWall::Ring { ring0, ring1 } => {
4272                    vec![if end == 0 {
4273                        ring0.clone()
4274                    } else {
4275                        ring1.clone()
4276                    }]
4277                }
4278                LoopWall::Chain { bottoms, tops } => {
4279                    if end == 0 {
4280                        bottoms.clone()
4281                    } else {
4282                        tops.clone()
4283                    }
4284                }
4285            });
4286        }
4287        let mut reach = 1.0_f64;
4288        for edges in &loop_edges {
4289            for edge in edges {
4290                let (curve, range) = spine_curve_of(model, edge)?;
4291                for t in 0..8 {
4292                    let p =
4293                        curve.point_at(range.0 + (range.1 - range.0) * f64::from(t) / 8.0, tol)?;
4294                    reach = reach.max(p.distance(at) * 2.0);
4295                }
4296            }
4297        }
4298        let cap_surface: SurfaceGeometry =
4299            PlaneSurface::over(cap_plane, (-reach, reach), (-reach, reach))?.into();
4300        let mut wires: Vec<Shape> = Vec::with_capacity(loop_edges.len());
4301        for edges in &loop_edges {
4302            wires.push(ogeom_algo::make_wire(model, edges, tol)?.shape);
4303        }
4304        let face = ogeom_algo::make_face(model, cap_surface.clone(), &wires, tol)?.shape;
4305        let cap_id = {
4306            let Some(node) = model.node(&face) else {
4307                ogeom_bail!(Dangling, "the cap just built is not in this model");
4308            };
4309            let ogeom_topo::NodeData::Face(data) = node.data() else {
4310                ogeom_bail!(Construction, "the cap holds no face data");
4311            };
4312            data.surface
4313        };
4314        let frame = cap_plane.frame();
4315        for edges in &loop_edges {
4316            for edge in edges {
4317                let (curve, range) = spine_curve_of(model, edge)?;
4318                let ogeom_geom::Curve::BSpline(bs) = &curve else {
4319                    ogeom_bail!(Construction, "a swept ring is not a spline");
4320                };
4321                // A planar polynomial spline's chart image is the same-degree
4322                // spline of the projected control points: affine, so exact.
4323                let control2: Vec<Point2> = bs
4324                    .control_points()
4325                    .iter()
4326                    .map(|w| {
4327                        let local = frame.to_local(w.point());
4328                        Point2::new(local.x, local.y)
4329                    })
4330                    .collect();
4331                let pcurve: ogeom_geom::PlanarCurve =
4332                    ogeom_geom::BSpline2d::new(bs.knots().clone(), control2, tol)?.into();
4333                ogeom_algo::attach_pcurve(
4334                    model,
4335                    edge,
4336                    pcurve,
4337                    cap_id,
4338                    ogeom_topo::Location::identity(),
4339                    range,
4340                )?;
4341            }
4342        }
4343        faces.push(face);
4344    }
4345
4346    let sewn = sew(model, &faces, tol)?;
4347    let mut built = if ring {
4348        // A holed ring sews into one shell per profile loop: the outer
4349        // bounds the material, each hole a void tunnel. Largest bound
4350        // first, the voids' faces already turned at build.
4351        if sewn.shells.is_empty() {
4352            ogeom_bail!(Construction, "the pipe shell did not close");
4353        }
4354        for shell in &sewn.shells {
4355            if !ogeom_algo::is_shell_closed(model, shell)? {
4356                ogeom_bail!(Construction, "the pipe shell did not close");
4357            }
4358        }
4359        let mut ordered = sewn.shells.clone();
4360        let mut sized: Vec<(f64, Shape)> = Vec::with_capacity(ordered.len());
4361        for shell in ordered.drain(..) {
4362            let bound = ogeom_algo::shape_bounds(model, &shell, tol)?;
4363            sized.push((bound.diagonal(), shell));
4364        }
4365        sized.sort_by(|a, b| b.0.partial_cmp(&a.0).unwrap_or(core::cmp::Ordering::Equal));
4366        let shells: Vec<Shape> = sized.into_iter().map(|(_, s)| s).collect();
4367        make_solid(model, &shells)?
4368    } else {
4369        if sewn.shells.len() != 1 || !ogeom_algo::is_shell_closed(model, &sewn.shells[0])? {
4370            if std::env::var_os("OGEOM_DEBUG_SWEEP").is_some() {
4371                eprintln!(
4372                    "SWEEP: {} shells, {} free edges from {} faces",
4373                    sewn.shells.len(),
4374                    sewn.free_edges.len(),
4375                    faces.len()
4376                );
4377                for edge in &sewn.free_edges {
4378                    let (curve, range) = spine_curve_of(model, edge)?;
4379                    let a = curve.point_at(range.0, tol)?;
4380                    let b = curve.point_at(range.1, tol)?;
4381                    let m = curve.point_at(f64::midpoint(range.0, range.1), tol)?;
4382                    let t = model.tolerance_of(edge)?.map_or(0.0, |t| t.get());
4383                    eprintln!(
4384                        "  free ({:.3},{:.3},{:.3}) -> ({:.3},{:.3},{:.3}) via ({:.3},{:.3},{:.3}) tol {t:.2e}",
4385                        a.x, a.y, a.z, b.x, b.y, b.z, m.x, m.y, m.z
4386                    );
4387                }
4388            }
4389            ogeom_bail!(Construction, "the pipe shell did not close");
4390        }
4391        make_solid(model, std::slice::from_ref(&sewn.shells[0]))?
4392    };
4393    built.history.generate(profile, built.shape.clone());
4394    built.history.generate(spine, built.shape.clone());
4395    Ok(built)
4396}
4397
4398/// An edge's 3D curve and range, cloned out of the model.
4399/// A pipe shell whose spine turns a skew corner against a curved leg,
4400/// built as pieces between such corners and fused.
4401///
4402/// With the frame reflected across the mitre plane, a straight leg's walls
4403/// and a curved leg's cut that plane in sections that differ on the inside
4404/// of the turn: no single join row closes both. Each piece is swept on
4405/// straight past its corners (the frame carries unchanged along a straight
4406/// run), trimmed by each corner's mitre plane, and the pieces fused: where
4407/// their sections on the plane coincide the caps melt, and where they
4408/// differ the difference stands as a face of the mitre plane. Exact, and
4409/// the plain mitre wherever the two sections agree.
4410#[allow(clippy::too_many_arguments, reason = "one construction, all its data")]
4411fn mitred_pieces(
4412    model: &mut Model,
4413    profile: &Shape,
4414    spine: &Shape,
4415    stations: &[SpineStation],
4416    skew: &[(usize, usize)],
4417    ring: bool,
4418    probes: &[(f64, f64)],
4419    tolerance: f64,
4420    tol: Tolerances,
4421) -> OgeomResult<Built> {
4422    if model.kind_of(profile)? == ShapeType::Wire {
4423        // A closed planar wire sweeps the walls of the face it bounds: the
4424        // solid pieces are mitred and fused as for that face, and its end
4425        // caps (on the planes square to the spine's ends) are taken off.
4426        if !ogeom_algo::is_wire_closed(model, profile, tol)? {
4427            ogeom_bail!(
4428                Construction,
4429                "a skew corner against a curved leg is mitred by fusing solid \
4430                 pieces; an open wire bounds no face to sweep round it"
4431            );
4432        }
4433        let Some(plane) = ogeom_algo::find_plane(model, profile, tol)? else {
4434            ogeom_bail!(Construction, "a pipe shell sweeps a planar profile");
4435        };
4436        let reach = 1e4_f64;
4437        let surface: SurfaceGeometry =
4438            PlaneSurface::over(plane, (-reach, reach), (-reach, reach))?.into();
4439        let face = ogeom_algo::make_face(model, surface, std::slice::from_ref(profile), tol)?.shape;
4440        let face = realized_profile(model, &face, &Transform::IDENTITY, tol)?;
4441        let solid = mitred_pieces(
4442            model, &face, spine, stations, skew, ring, probes, tolerance, tol,
4443        )?
4444        .shape;
4445        let ends: Vec<(Point, Vector)> = if ring {
4446            Vec::new()
4447        } else {
4448            vec![
4449                (stations[0].at, stations[0].tangent),
4450                (
4451                    stations[stations.len() - 1].at,
4452                    stations[stations.len() - 1].tangent,
4453                ),
4454            ]
4455        };
4456        let mut walls = Vec::new();
4457        for f in explore(model, &solid, Filter::OfType(ShapeType::Face))? {
4458            let Some(ogeom_topo::NodeData::Face(data)) = model.node(&f).map(|n| n.data()) else {
4459                continue;
4460            };
4461            let cap = match model.geometry().surface(data.surface) {
4462                Some(SurfaceGeometry::Plane(p)) => {
4463                    let placed = p.plane();
4464                    ends.iter().any(|(at, n)| {
4465                        placed.normal().vector().cross(*n).magnitude() <= tol.angular()
4466                            && placed.distance_to(*at) <= tol.confusion() * 100.0
4467                    })
4468                }
4469                _ => false,
4470            };
4471            if !cap {
4472                walls.push(f);
4473            }
4474        }
4475        let sewn = sew(model, &walls, tol)?;
4476        let shape = match sewn.shells.as_slice() {
4477            [shell] => shell.clone(),
4478            _ => ogeom_algo::make_compound(model, &sewn.shells)?.shape,
4479        };
4480        let mut history = History::new();
4481        history.generate(profile, shape.clone());
4482        history.generate(spine, shape.clone());
4483        return Ok(Built::new(shape, history));
4484    }
4485    if model.kind_of(profile)? != ShapeType::Face {
4486        ogeom_bail!(Construction, "a pipe shell sweeps a planar wire or face");
4487    }
4488    let edges: Vec<Shape> = match model.kind_of(spine)? {
4489        ShapeType::Edge => vec![spine.clone()],
4490        _ => model.ordered_children_of(spine)?,
4491    };
4492    let normals = rmf_normals(stations);
4493    let reach_out = probes
4494        .iter()
4495        .map(|(a, b)| a.hypot(*b))
4496        .fold(0.0_f64, f64::max);
4497
4498    // Each split corner: the junction's edges, its point, the tangents
4499    // either side, the frame the far side starts in, and how far each side
4500    // runs on past it to cover the mitre plane across the whole profile.
4501    struct Split {
4502        edge_before: usize,
4503        at: Point,
4504        before: Vector,
4505        after: Vector,
4506        frame_after: Vector,
4507        run_on: f64,
4508    }
4509    let mut splits: Vec<Split> = Vec::with_capacity(skew.len());
4510    for &(k, next) in skew {
4511        let (before, after) = (stations[k].tangent, stations[next].tangent);
4512        let turn = before.dot(after).clamp(-1.0, 1.0).acos();
4513        let half = (turn * 0.5).cos();
4514        if half < 0.05 {
4515            ogeom_bail!(
4516                Construction,
4517                "the spine all but doubles back at a corner; a mitre there \
4518                 runs off to infinity"
4519            );
4520        }
4521        splits.push(Split {
4522            edge_before: stations[k].edge,
4523            at: stations[k].at,
4524            before,
4525            after,
4526            frame_after: normals[next],
4527            // The mitre plane stands at most `R·tan(φ/2)` past the corner
4528            // along either leg for a profile reaching `R` from the spine;
4529            // half as far again clears it with room.
4530            run_on: reach_out * ((turn * 0.5).tan() * 1.5 + 0.1),
4531        });
4532    }
4533    splits.sort_by_key(|s| s.edge_before);
4534
4535    // Pieces as runs of spine edges, each between split corners (or the
4536    // open spine's own ends).
4537    let count = edges.len();
4538    let mut pieces: Vec<(Vec<usize>, Option<usize>, Option<usize>)> = Vec::new();
4539    if ring {
4540        for (i, split) in splits.iter().enumerate() {
4541            let next = &splits[(i + 1) % splits.len()];
4542            let mut run = Vec::new();
4543            let mut e = (split.edge_before + 1) % count;
4544            loop {
4545                run.push(e);
4546                if e == next.edge_before {
4547                    break;
4548                }
4549                e = (e + 1) % count;
4550            }
4551            pieces.push((run, Some(i), Some((i + 1) % splits.len())));
4552        }
4553    } else {
4554        let mut first = 0;
4555        for (i, split) in splits.iter().enumerate() {
4556            pieces.push((
4557                (first..=split.edge_before).collect(),
4558                i.checked_sub(1),
4559                Some(i),
4560            ));
4561            first = split.edge_before + 1;
4562        }
4563        pieces.push(((first..count).collect(), splits.len().checked_sub(1), None));
4564    }
4565
4566    let x0 = normals[0];
4567    let start_frame = Frame::new(
4568        stations[0].at,
4569        Direction::new(stations[0].tangent, tol)?,
4570        Direction::new(x0, tol)?,
4571        tol,
4572    )?;
4573    // One block per split, standing on the mitre plane on the far side of
4574    // the corner: the piece before the corner is cut by it and the piece
4575    // after keeps what it shares with it, so both sides' caps are pieces
4576    // of the block's one face, on one surface and one chart, which is what
4577    // lets the fuse melt them.
4578    let mut blocks: Vec<Shape> = Vec::with_capacity(splits.len());
4579    for split in &splits {
4580        let n = (split.before + split.after) / (split.before + split.after).magnitude();
4581        let normal = Direction::new(n, tol)?;
4582        let plane = Plane::through(split.at, normal);
4583        let reach = (split.run_on + reach_out) * 4.0;
4584        let frame = plane.frame();
4585        let (u, v) = (frame.x().vector(), frame.y().vector());
4586        let corners: Vec<Point> = [(-1.0, -1.0), (1.0, -1.0), (1.0, 1.0), (-1.0, 1.0)]
4587            .iter()
4588            .map(|(a, b)| split.at + u * (a * reach) + v * (b * reach))
4589            .collect();
4590        let wire = ogeom_algo::make_polygon(model, &corners, true, tol)?.shape;
4591        let edges = explore(model, &wire, Filter::OfType(ShapeType::Edge))?;
4592        let surface: SurfaceGeometry = PlaneSurface::over(
4593            plane,
4594            (-reach * 2.0, reach * 2.0),
4595            (-reach * 2.0, reach * 2.0),
4596        )?
4597        .into();
4598        let base = ogeom_algo::make_face_with_pcurves(model, surface, &[edges], tol)?.shape;
4599        let block = ogeom_algo::make_prism(model, &base, n * (reach * 2.0), tol)?.shape;
4600        blocks.push(block);
4601    }
4602
4603    let mut result: Option<Shape> = None;
4604    for (run, start, end) in pieces {
4605        let mut wire_edges: Vec<Shape> = Vec::new();
4606        let traversal = |model: &Model, e: usize, at_start: bool| -> OgeomResult<Shape> {
4607            // In the spine's own sense: a reversed edge's vertices come
4608            // back end first.
4609            let Some((a, b)) = ogeom_algo::edge_vertices(model, &edges[e])? else {
4610                ogeom_bail!(Construction, "a spine edge has no vertices");
4611            };
4612            Ok(if at_start { a } else { b })
4613        };
4614        if let Some(i) = start {
4615            let split = &splits[i];
4616            let far = ogeom_algo::make_vertex(model, split.at - split.after * split.run_on).shape;
4617            let near = traversal(model, run[0], true)?;
4618            let line: ogeom_geom::Curve =
4619                LineCurve::segment(split.at - split.after * split.run_on, split.at, tol)?.into();
4620            let domain = line.domain();
4621            wire_edges
4622                .push(ogeom_algo::make_edge_between(model, line, domain, &far, &near, tol)?.shape);
4623        }
4624        wire_edges.extend(run.iter().map(|&e| edges[e].clone()));
4625        if let Some(i) = end {
4626            let split = &splits[i];
4627            let near = traversal(model, run[run.len() - 1], false)?;
4628            let far = ogeom_algo::make_vertex(model, split.at + split.before * split.run_on).shape;
4629            let line: ogeom_geom::Curve =
4630                LineCurve::segment(split.at, split.at + split.before * split.run_on, tol)?.into();
4631            let domain = line.domain();
4632            wire_edges
4633                .push(ogeom_algo::make_edge_between(model, line, domain, &near, &far, tol)?.shape);
4634        }
4635        let sub_spine = ogeom_algo::make_wire(model, &wire_edges, tol)?.shape;
4636        // The profile where this piece starts: the spine's own start keeps
4637        // the caller's; a piece starting past a corner takes the profile
4638        // moved into the frame the corner's far side starts in, set back
4639        // along its run-on.
4640        let placed = match start {
4641            None => profile.clone(),
4642            Some(i) => {
4643                let split = &splits[i];
4644                let target = Frame::new(
4645                    split.at - split.after * split.run_on,
4646                    Direction::new(split.after, tol)?,
4647                    Direction::new(split.frame_after, tol)?,
4648                    tol,
4649                )?;
4650                let motion = Transform::from_frame(&target) * Transform::to_frame(&start_frame);
4651                realized_profile(model, profile, &motion, tol)?
4652            }
4653        };
4654        let mut piece = make_pipe_shell(model, &placed, &sub_spine, false, tolerance, tol)?.shape;
4655        if model.kind_of(&piece)? != ShapeType::Solid {
4656            ogeom_bail!(Construction, "a mitred piece did not sweep into a solid");
4657        }
4658        if let Some(i) = start {
4659            piece = ogeom_bool::common(model, &piece, &blocks[i], tol)?.shape;
4660        }
4661        if let Some(i) = end {
4662            piece = ogeom_bool::cut(model, &piece, &blocks[i], tol)?.shape;
4663        }
4664        result = Some(match result {
4665            None => piece,
4666            Some(held) => ogeom_bool::fuse(model, &held, &piece, tol)?.shape,
4667        });
4668    }
4669    let Some(shape) = result else {
4670        ogeom_bail!(Construction, "the spine produced no piece to sweep");
4671    };
4672    let mut history = History::new();
4673    history.generate(profile, shape.clone());
4674    for edge in &edges {
4675        history.generate(edge, shape.clone());
4676    }
4677    Ok(Built::new(shape, history))
4678}
4679
4680/// A planar profile face rebuilt under a rigid motion: every edge's curve
4681/// moved and re-bounded, vertices shared, the face on the moved plane.
4682fn realized_profile(
4683    model: &mut Model,
4684    profile: &Shape,
4685    motion: &Transform,
4686    tol: Tolerances,
4687) -> OgeomResult<Shape> {
4688    realized_profile_wound(model, profile, motion, None, tol)
4689}
4690
4691/// As [`realized_profile`], each ring wound about `about` where given: the
4692/// outer ring turning positively, every hole the other way.
4693fn realized_profile_wound(
4694    model: &mut Model,
4695    profile: &Shape,
4696    motion: &Transform,
4697    about: Option<Vector>,
4698    tol: Tolerances,
4699) -> OgeomResult<Shape> {
4700    use ogeom_geom::Transformable as _;
4701    let Some(plane) = ogeom_algo::find_plane(model, profile, tol)? else {
4702        ogeom_bail!(Construction, "a pipe shell sweeps a planar profile");
4703    };
4704    let moved_plane = Plane::through(
4705        motion.apply(plane.origin()),
4706        Direction::new(motion.apply_vector(plane.normal().vector()), tol)?,
4707    );
4708    let mut vertices: std::collections::HashMap<ogeom_topo::TShapeId, Shape> =
4709        std::collections::HashMap::new();
4710    let mut edge_copies: std::collections::HashMap<ogeom_topo::TShapeId, Shape> =
4711        std::collections::HashMap::new();
4712    let mut wires: Vec<Vec<Shape>> = Vec::new();
4713    for wire in explore(model, profile, Filter::OfType(ShapeType::Wire))? {
4714        let mut ring = Vec::new();
4715        for edge in model.ordered_children_of(&wire)? {
4716            let copy = match edge_copies.get(&edge.node()) {
4717                Some(done) => done.clone(),
4718                None => {
4719                    let (curve, range) = spine_curve_of(model, &edge)?;
4720                    let placed = curve.transformed(&edge.transform(model.datums())?, tol)?;
4721                    let moved = placed.transformed(motion, tol)?;
4722                    // The copy is of the edge itself, in its own sense; the
4723                    // ring's use of it is reapplied below.
4724                    let own = if edge.orientation() == ogeom_topo::Orientation::Reversed {
4725                        edge.reversed()
4726                    } else {
4727                        edge.clone()
4728                    };
4729                    let Some((a, b)) = ogeom_algo::edge_vertices(model, &own)? else {
4730                        ogeom_bail!(Construction, "a profile edge has no vertices");
4731                    };
4732                    let mut ends = Vec::with_capacity(2);
4733                    for v in [a, b] {
4734                        let key = v.node();
4735                        let held = match vertices.get(&key) {
4736                            Some(done) => done.clone(),
4737                            None => {
4738                                let Some(data) = model.node(&v).and_then(|n| n.data().as_vertex())
4739                                else {
4740                                    ogeom_bail!(Construction, "a profile vertex holds no data");
4741                                };
4742                                let at = v.transform(model.datums())?.apply(data.point);
4743                                let fresh = ogeom_algo::make_vertex(model, motion.apply(at)).shape;
4744                                vertices.insert(key, fresh.clone());
4745                                fresh
4746                            }
4747                        };
4748                        ends.push(held);
4749                    }
4750                    let fresh = ogeom_algo::make_edge_between(
4751                        model, moved, range, &ends[0], &ends[1], tol,
4752                    )?
4753                    .shape;
4754                    edge_copies.insert(edge.node(), fresh.clone());
4755                    fresh
4756                }
4757            };
4758            ring.push(if edge.orientation() == ogeom_topo::Orientation::Reversed {
4759                copy.reversed()
4760            } else {
4761                copy
4762            });
4763        }
4764        if let Some(axis) = about {
4765            let turning = ring_turning(model, &ring, axis, tol)?;
4766            let outer = wires.is_empty();
4767            if (turning > 0.0) != outer {
4768                ring = ring.iter().rev().map(Shape::reversed).collect();
4769            }
4770        }
4771        wires.push(ring);
4772    }
4773    let reach = 1e4_f64;
4774    let surface: SurfaceGeometry =
4775        PlaneSurface::over(moved_plane, (-reach, reach), (-reach, reach))?.into();
4776    Ok(ogeom_algo::make_face_with_pcurves(model, surface, &wires, tol)?.shape)
4777}
4778
4779/// Twice the signed area a ring of edges encloses about `axis`, from its
4780/// edges sampled in the ring's own sense.
4781fn ring_turning(model: &Model, ring: &[Shape], axis: Vector, tol: Tolerances) -> OgeomResult<f64> {
4782    let mut points: Vec<Point> = Vec::new();
4783    for edge in ring {
4784        let (curve, range) = spine_curve_of(model, edge)?;
4785        let reversed = edge.orientation() == ogeom_topo::Orientation::Reversed;
4786        for i in 0..32 {
4787            let f = f64::from(i) / 32.0;
4788            let t = if reversed {
4789                range.1 - (range.1 - range.0) * f
4790            } else {
4791                range.0 + (range.1 - range.0) * f
4792            };
4793            points.push(curve.point_at(t, tol)?);
4794        }
4795    }
4796    let Some(&origin) = points.first() else {
4797        return Ok(0.0);
4798    };
4799    let n = points.len();
4800    Ok((0..n)
4801        .map(|i| {
4802            (points[i] - origin)
4803                .cross(points[(i + 1) % n] - origin)
4804                .dot(axis)
4805        })
4806        .sum())
4807}
4808
4809fn spine_curve_of(model: &Model, edge: &Shape) -> OgeomResult<(ogeom_geom::Curve, (f64, f64))> {
4810    let Some(data) = model.node(edge).and_then(|n| n.data().as_edge()) else {
4811        ogeom_bail!(Construction, "an edge holds no data");
4812    };
4813    let Some(EdgeRepr::Curve3d { curve, range, .. }) = data.curve3d() else {
4814        ogeom_bail!(Construction, "an edge has no curve");
4815    };
4816    let Some(geometry) = model.geometry().curve(*curve) else {
4817        ogeom_bail!(Dangling, "curve is not in this model");
4818    };
4819    Ok((geometry.clone(), *range))
4820}
4821
4822/// The pipe shell around a spine that loops back on itself: one wall,
4823/// closed both ways round, no caps at all.
4824///
4825/// The frames are rotation-minimizing with the loop's holonomy paid off:
4826/// transported round a closed spine, the frame comes home twisted by some
4827/// angle, and that twist is spread back along the arc so the last station's
4828/// frame *is* the first's; without it the closed fit fights a helical
4829/// grid. The profile may be smooth or faceted, and each hole sweeps a void
4830/// tunnel of its own shell; the Frenet law rides the loop as well.
4831fn closed_pipe_shell(
4832    model: &mut Model,
4833    profile: &Shape,
4834    spine: &Shape,
4835    mut stations: Vec<SpineStation>,
4836    frenet: bool,
4837    tolerance: f64,
4838    tol: Tolerances,
4839) -> OgeomResult<Built> {
4840    // The walk visits the join twice; the loop owns it once.
4841    stations.pop();
4842    if stations.len() < 3 {
4843        ogeom_bail!(Construction, "a closed spine needs room to turn");
4844    }
4845    // A sharp corner is a kink, and a kinked ring is mitred by the caller;
4846    // one arriving here has a heading that jumps between two stations.
4847    for i in 0..stations.len() {
4848        let next = &stations[(i + 1) % stations.len()];
4849        if stations[i].tangent.dot(next.tangent) < 0.9 {
4850            ogeom_bail!(
4851                Construction,
4852                "a closed spine turns too sharply between two of its stations to skin"
4853            );
4854        }
4855    }
4856
4857    let loops: Vec<Shape> = match model.kind_of(profile)? {
4858        ShapeType::Face => explore(model, profile, Filter::OfType(ShapeType::Wire))?,
4859        ShapeType::Wire => vec![profile.clone()],
4860        other => ogeom_bail!(
4861            Construction,
4862            "a pipe shell sweeps a planar wire or face, not a {other:?}"
4863        ),
4864    };
4865    // Every wire sweeps its own closed shell: the outer boundary first,
4866    // each hole a void tunnel inside it.
4867    let profile_loop = &loops[0];
4868    let edges = model.ordered_children_of(profile_loop)?;
4869    let smooth = edges.len() == 1
4870        && ogeom_algo::edge_vertices(model, &edges[0])?.is_some_and(|(a, b)| a.is_same(&b));
4871    let Some(plane) = ogeom_algo::find_plane(model, profile, tol)? else {
4872        ogeom_bail!(Construction, "a pipe shell sweeps a planar profile");
4873    };
4874    let t0 = stations[0].tangent;
4875    // Square to the spine's own start tangent, read exactly off its first
4876    // edge, to an angle's tolerance.
4877    let exact_t0 = {
4878        let first = match model.kind_of(spine)? {
4879            ShapeType::Edge => spine.clone(),
4880            _ => model.ordered_children_of(spine)?[0].clone(),
4881        };
4882        let (curve, range) = spine_curve_of(model, &first)?;
4883        let reversed = first.orientation() == ogeom_topo::Orientation::Reversed;
4884        let d = curve.d1_at(if reversed { range.1 } else { range.0 }, tol)?;
4885        let d = if reversed { -d } else { d };
4886        d / d.magnitude()
4887    };
4888    if plane.normal().vector().cross(exact_t0).magnitude() > tol.angular() {
4889        ogeom_bail!(
4890            Construction,
4891            "the profile leans along its spine; a pipe shell runs square to \
4892             the start"
4893        );
4894    }
4895    if plane.distance_to(stations[0].at) > tol.confusion() * 100.0 {
4896        ogeom_bail!(
4897            Construction,
4898            "the profile does not sit at the spine's start"
4899        );
4900    }
4901
4902    // Frames with the loop's mismatch paid off: carry once more back to the
4903    // start, read the twist between departure and return, and spread it
4904    // along the arc. Rotation-minimizing frames owe this for their
4905    // holonomy; the Frenet law owes it too, because straight stretches
4906    // carry the frame through by continuation and the continuation is
4907    // path-dependent. One reconciliation serves both.
4908    let mut normals: Vec<Vector> = if frenet {
4909        // The Frenet frame is the spine's own, single-valued round a loop:
4910        // read with wrapped neighbours it closes on itself and owes no
4911        // reconciliation. Read from a walk that visits the join twice it
4912        // does not: the one-sided differences at the walk's two ends
4913        // disagree with the interior, and the strips built on them miss
4914        // each other at the join by that kink.
4915        frenet_normals_closed(&stations, tol)?
4916    } else {
4917        let mut extended = stations.clone();
4918        extended.push(stations[0]);
4919        let carried = rmf_normals(&extended);
4920        let (n0, n_home) = (carried[0], carried[carried.len() - 1]);
4921        let twist = (n0.cross(n_home).dot(t0)).atan2(n0.dot(n_home));
4922        let mut lengths = vec![0.0_f64];
4923        for pair in extended.windows(2) {
4924            let last = lengths[lengths.len() - 1];
4925            lengths.push(last + pair[0].at.distance(pair[1].at));
4926        }
4927        let total = lengths[lengths.len() - 1];
4928        carried
4929            .iter()
4930            .take(stations.len())
4931            .enumerate()
4932            .map(|(i, n)| {
4933                let phi = -twist * lengths[i] / total;
4934                let t = extended[i].tangent;
4935                *n * phi.cos() + t.cross(*n) * phi.sin()
4936            })
4937            .collect()
4938    };
4939    for (n, station) in normals.iter_mut().zip(&stations) {
4940        // Re-square each corrected normal against its own tangent.
4941        let v = *n - station.tangent * n.dot(station.tangent);
4942        *n = v / v.magnitude();
4943    }
4944
4945    let x0 = normals[0];
4946    let origin = stations[0].at;
4947    let mut shells: Vec<Shape> = Vec::with_capacity(loops.len());
4948    for (li, wire) in loops.iter().enumerate() {
4949        let wire_edges = model.ordered_children_of(wire)?;
4950        let wire_smooth = wire_edges.len() == 1
4951            && ogeom_algo::edge_vertices(model, &wire_edges[0])?
4952                .is_some_and(|(a, b)| a.is_same(&b));
4953        let shell = closed_loop_shell(
4954            model,
4955            wire,
4956            &wire_edges,
4957            wire_smooth,
4958            &stations,
4959            &normals,
4960            (origin, x0),
4961            tolerance,
4962            tol,
4963        )?;
4964        // The material side follows each loop's own winding, so it is
4965        // read off the shell itself: the outer shell faces out of what it
4966        // encloses, a void's faces toward its tunnel.
4967        let enclosed = shell_signed_volume(model, &shell, tol)?;
4968        let outward = enclosed > 0.0;
4969        shells.push(if outward == (li == 0) {
4970            shell
4971        } else {
4972            shell.reversed()
4973        });
4974    }
4975    let mut built = make_solid(model, &shells)?;
4976    built.history.generate(profile, built.shape.clone());
4977    built.history.generate(spine, built.shape.clone());
4978    let _ = (smooth, profile_loop, edges);
4979    Ok(built)
4980}
4981
4982/// The volume a closed shell encloses as it faces, from its mesh: negative
4983/// where its faces point into what it bounds.
4984fn shell_signed_volume(model: &Model, shell: &Shape, tol: Tolerances) -> OgeomResult<f64> {
4985    let mesh = ogeom_mesh::triangulate(model, shell, ogeom_mesh::Deflection::default(), tol)?;
4986    Ok(mesh
4987        .triangles
4988        .iter()
4989        .map(|t| {
4990            let [a, b, c] = t.map(|i| mesh.positions[i as usize].to_vector());
4991            a.dot(b.cross(c)) / 6.0
4992        })
4993        .sum())
4994}
4995
4996/// Sample a spine (one edge or a wire of them) into stations, each edge
4997/// given a station count by its own turning.
4998fn shell_stations(model: &Model, spine: &Shape, tol: Tolerances) -> OgeomResult<Vec<SpineStation>> {
4999    let edges: Vec<Shape> = match model.kind_of(spine)? {
5000        ShapeType::Edge => vec![spine.clone()],
5001        ShapeType::Wire => model.ordered_children_of(spine)?,
5002        other => ogeom_bail!(
5003            Construction,
5004            "a pipe shell runs along an edge or a wire, not a {other:?}"
5005        ),
5006    };
5007    if edges.is_empty() {
5008        ogeom_bail!(Construction, "the spine has no edge to run along");
5009    }
5010    let mut stations: Vec<SpineStation> = Vec::new();
5011    for (ei, edge) in edges.iter().enumerate() {
5012        let (curve, range) = {
5013            let Some(data) = model.node(edge).and_then(|n| n.data().as_edge()) else {
5014                ogeom_bail!(Construction, "a spine edge holds no data");
5015            };
5016            let Some(EdgeRepr::Curve3d { curve, range, .. }) = data.curve3d() else {
5017                ogeom_bail!(Construction, "a spine edge has no curve");
5018            };
5019            let Some(geometry) = model.geometry().curve(*curve) else {
5020                ogeom_bail!(Dangling, "curve is not in this model");
5021            };
5022            (geometry.clone(), *range)
5023        };
5024        let reversed = edge.orientation() == ogeom_topo::Orientation::Reversed;
5025        // Stations by turning: sample tangents coarsely, sum the angles, and
5026        // give each edge enough stations that no step turns more than a few
5027        // degrees. A straight edge keeps a healthy minimum for the fit.
5028        let turning = {
5029            let mut sum = 0.0_f64;
5030            let mut last: Option<Vector> = None;
5031            for i in 0..=16 {
5032                let t = range.0 + (range.1 - range.0) * f64::from(i) / 16.0;
5033                let d = curve.d1_at(t, tol)?;
5034                let m = d.magnitude();
5035                if m <= tol.confusion() {
5036                    continue;
5037                }
5038                let u = d / m;
5039                if let Some(prev) = last {
5040                    sum += prev.dot(u).clamp(-1.0, 1.0).acos() * 16.0 / 16.0;
5041                }
5042                last = Some(u);
5043            }
5044            sum
5045        };
5046        #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
5047        let count = (turning / (core::f64::consts::TAU / 64.0)).ceil().max(8.0) as usize;
5048        for i in 0..=count {
5049            #[allow(clippy::cast_precision_loss)]
5050            let f = (i as f64) / (count as f64);
5051            let t = if reversed {
5052                range.1 - (range.1 - range.0) * f
5053            } else {
5054                range.0 + (range.1 - range.0) * f
5055            };
5056            let p = curve.point_at(t, tol)?;
5057            let d = curve.d1_at(t, tol)?;
5058            let m = d.magnitude();
5059            if m <= tol.confusion() {
5060                ogeom_bail!(Construction, "the spine is degenerate at {t}");
5061            }
5062            let tangent = if reversed { -(d / m) } else { d / m };
5063            if let Some(prev) = stations.last()
5064                && prev.edge == ei
5065                && prev.at.distance(p) <= tol.confusion()
5066                && prev.tangent.cross(tangent).magnitude() <= tol.angular()
5067                && prev.tangent.dot(tangent) > 0.0
5068            {
5069                continue;
5070            }
5071            // A station coincident with the last but on the next edge is a
5072            // twin: heading elsewhere it is a *corner* the sweep mitres, and
5073            // heading on it is a smooth junction where the next edge's own
5074            // run of skin begins, since one fit across two curves' joins
5075            // cannot follow the step in their curvature.
5076            stations.push(SpineStation {
5077                at: p,
5078                tangent,
5079                edge: ei,
5080                t,
5081            });
5082        }
5083    }
5084    if stations.len() < 2 {
5085        ogeom_bail!(Construction, "the spine collapses to a point");
5086    }
5087    Ok(stations)
5088}
5089
5090/// Frenet normals: each station's frame turns with the spine's own
5091/// curvature, read from the tangents' finite differences. Straight runs
5092/// carry the last bending station's normal forward; a spine that never
5093/// bends has no Frenet frame at all and is refused by name.
5094fn frenet_normals(stations: &[SpineStation], tol: Tolerances) -> OgeomResult<Vec<Vector>> {
5095    let mut normals: Vec<Option<Vector>> = Vec::with_capacity(stations.len());
5096    for i in 0..stations.len() {
5097        let (before, after) = (
5098            &stations[i.saturating_sub(1)],
5099            &stations[(i + 1).min(stations.len() - 1)],
5100        );
5101        let dt = after.tangent - before.tangent;
5102        let t = stations[i].tangent;
5103        let bend = dt - t * dt.dot(t);
5104        let m = bend.magnitude();
5105        normals.push(if m > tol.angular().max(1e-9) {
5106            Some(bend / m)
5107        } else {
5108            None
5109        });
5110    }
5111    // Carry forward, then backward, so straight lead-ins take the first
5112    // bend's frame rather than none.
5113    let mut carried: Vec<Vector> = Vec::with_capacity(stations.len());
5114    let mut last: Option<Vector> = None;
5115    for n in &normals {
5116        if let Some(n) = n {
5117            last = Some(*n);
5118        }
5119        carried.push(last.unwrap_or(Vector::new(0.0, 0.0, 0.0)));
5120    }
5121    let mut ahead: Option<Vector> = None;
5122    for i in (0..stations.len()).rev() {
5123        if let Some(n) = normals[i] {
5124            ahead = Some(n);
5125        } else if carried[i].magnitude() < 0.5
5126            && let Some(n) = ahead
5127        {
5128            carried[i] = n;
5129        }
5130    }
5131    if carried.iter().any(|n| n.magnitude() < 0.5) {
5132        ogeom_bail!(
5133            Construction,
5134            "a straight spine has no Frenet frame; use the \
5135             rotation-minimizing default"
5136        );
5137    }
5138    Ok(carried)
5139}
5140
5141/// Frenet normals round a closed loop: each station's bend read from its
5142/// neighbours across the join as well, so the field is periodic. A loop
5143/// with a straight stretch carries the last bend's normal through it, as
5144/// the open form does; a loop that never bends has no Frenet frame.
5145fn frenet_normals_closed(stations: &[SpineStation], tol: Tolerances) -> OgeomResult<Vec<Vector>> {
5146    let n = stations.len();
5147    let mut normals: Vec<Option<Vector>> = Vec::with_capacity(n);
5148    for i in 0..n {
5149        let (before, after) = (&stations[(i + n - 1) % n], &stations[(i + 1) % n]);
5150        let dt = after.tangent - before.tangent;
5151        let t = stations[i].tangent;
5152        let bend = dt - t * dt.dot(t);
5153        let m = bend.magnitude();
5154        normals.push(if m > tol.angular().max(1e-9) {
5155            Some(bend / m)
5156        } else {
5157            None
5158        });
5159    }
5160    let Some(first_bend) = normals.iter().position(Option::is_some) else {
5161        ogeom_bail!(
5162            Construction,
5163            "a straight spine has no Frenet frame; use the \
5164             rotation-minimizing default"
5165        );
5166    };
5167    // Carry forward round the loop from the first bend, so a straight
5168    // stretch anywhere takes the bend behind it.
5169    let mut carried: Vec<Vector> = vec![Vector::new(0.0, 0.0, 0.0); n];
5170    let mut last = normals[first_bend].unwrap_or_else(|| unreachable!());
5171    for k in 0..n {
5172        let i = (first_bend + k) % n;
5173        if let Some(bend) = normals[i] {
5174            last = bend;
5175        }
5176        carried[i] = last;
5177    }
5178    Ok(carried)
5179}
5180
5181// --- the evolved shape -------------------------------------------------------
5182
5183/// One station of the spine: where it is, which way it runs, and how it gets
5184/// there.
5185struct Station {
5186    /// Where the traversal enters and leaves this edge.
5187    from: Point,
5188    to: Point,
5189    /// The unit tangent at each end, in the direction of travel.
5190    tangent_in: Vector,
5191    tangent_out: Vector,
5192    /// A straight run, or a turn about an axis through an angle.
5193    turn: Option<(ogeom_math::Axis, f64)>,
5194}
5195
5196/// Sweep a profile along a spine, the way a moulding runs round a frame.
5197///
5198/// The spine is a **planar** wire, or a planar face whose outer wire is taken.
5199/// The profile is a wire standing in a plane that contains the spine's own
5200/// normal, positioned where the spine starts. What comes back is what the
5201/// profile sweeps out as it travels the spine, always square to it:
5202///
5203/// - a straight spine edge extrudes the profile (a prism);
5204/// - a circular one turns it about that arc's own axis (a revolution);
5205/// - and each corner between them turns it about the corner, through exactly
5206///   the angle the spine turns there, which is the join the 2D offset makes
5207///   for the same reason.
5208///
5209/// Every piece is exact: the surfaces are the ones a prism and a revolution
5210/// give for the profile's own curves, and nothing is fitted. The pieces are
5211/// then unioned, which is the assembly's real name: consecutive pieces meet
5212/// on the *same* placed profile, and a coincident face is what the boolean
5213/// identifies rather than probes across.
5214///
5215/// # Volume or shell
5216///
5217/// The result is always a volume, and which spine is given is what says
5218/// whether there is one to have. A **closed** profile bounds its own section
5219/// and sweeps a solid along either kind of spine. An **open** one does not,
5220/// and there is exactly one honest way to close it: against the plane a
5221/// **face** spine was drawn in, whose own plane the profile's two ends must
5222/// reach. An open profile along a wire spine is refused, and the refusal says
5223/// which spine would close it.
5224///
5225/// # Errors
5226///
5227/// [`OgeomError::Construction`](ogeom_core::OgeomError::Construction) if the spine is
5228/// not a planar wire or face, if it carries an edge that is neither straight
5229/// nor circular, if the profile is not planar, if the profile's plane does not
5230/// contain the spine's normal or does not cut across it (a profile that leans
5231/// or lies along is not square to the spine), or if an open profile has no
5232/// spine plane to close against.
5233/// [`OgeomError::NotDone`](ogeom_core::OgeomError::NotDone) where a corner's turn
5234/// would sweep the profile across the corner itself, which no revolution can
5235/// express.
5236pub fn make_evolved(
5237    model: &mut Model,
5238    spine: &Shape,
5239    profile: &Shape,
5240    tol: Tolerances,
5241) -> OgeomResult<Built> {
5242    use ogeom_algo::{make_prism, make_revolution, transformed};
5243
5244    let (wire, capped_by_spine_plane) = match model.kind_of(spine)? {
5245        ShapeType::Face => {
5246            let wires = explore(model, spine, Filter::OfType(ShapeType::Wire))?;
5247            let Some(outer) = wires.first().cloned() else {
5248                ogeom_bail!(Construction, "a face with no wire has no spine to run");
5249            };
5250            (outer, true)
5251        }
5252        ShapeType::Wire => (spine.clone(), false),
5253        other => ogeom_bail!(
5254            Construction,
5255            "a {other:?} is not a spine; sweep along a wire or a planar face"
5256        ),
5257    };
5258
5259    let stations = spine_stations(model, &wire, tol)?;
5260    if stations.is_empty() {
5261        ogeom_bail!(Construction, "a spine with no edges goes nowhere");
5262    }
5263    let normal = spine_normal(&stations, tol)?;
5264    let (profile_origin, profile_normal) = profile_plane(model, profile, tol)?;
5265    if profile_normal.dot(normal.vector()).abs() > tol.angular() {
5266        ogeom_bail!(
5267            Construction,
5268            "the profile's plane must contain the spine's normal, or the \
5269             profile is not square to the spine it travels"
5270        );
5271    }
5272    let start = &stations[0];
5273    if profile_normal.cross(start.tangent_in).magnitude() > tol.angular() {
5274        ogeom_bail!(
5275            Construction,
5276            "the profile's plane must cut the spine across, not run along it: \
5277             the profile is not square to the spine it travels"
5278        );
5279    }
5280    let _ = profile_origin;
5281    let reference = (start.from, start.tangent_in);
5282
5283    // The profile as a face, which is what makes each swept piece a *solid*
5284    // and the assembly a union rather than a hopeful sew. An open profile is
5285    // closed against the spine's own plane, which is exactly what a face
5286    // spine offers and a wire spine does not.
5287    let section = profile_face(
5288        model,
5289        profile,
5290        profile_normal,
5291        capped_by_spine_plane.then(|| Plane::through(start.from, normal)),
5292        tol,
5293    )?;
5294
5295    let mut pieces: Vec<Shape> = Vec::new();
5296    for (index, station) in stations.iter().enumerate() {
5297        // The corner *before* this station, so the pieces come out in the
5298        // order the spine runs them.
5299        if index > 0 {
5300            let previous = &stations[index - 1];
5301            if let Some(piece) = corner_piece(
5302                model,
5303                &section,
5304                reference,
5305                previous.to,
5306                previous.tangent_out,
5307                station.tangent_in,
5308                normal,
5309                tol,
5310            )? {
5311                pieces.push(piece);
5312            }
5313        }
5314        let placed = transformed(
5315            model,
5316            &section,
5317            station_transform(reference, station.from, station.tangent_in, normal, tol)?,
5318        )?
5319        .shape;
5320        pieces.push(match station.turn {
5321            None => make_prism(model, &placed, station.to - station.from, tol)?.shape,
5322            Some((axis, angle)) => make_revolution(model, &placed, axis, angle, tol)?.shape,
5323        });
5324    }
5325    // A closed spine turns at the join between its last edge and its first
5326    // just as it does anywhere else.
5327    let last = &stations[stations.len() - 1];
5328    if last.to.distance(start.from) <= tol.confusion()
5329        && let Some(piece) = corner_piece(
5330            model,
5331            &section,
5332            reference,
5333            last.to,
5334            last.tangent_out,
5335            start.tangent_in,
5336            normal,
5337            tol,
5338        )?
5339    {
5340        pieces.push(piece);
5341    }
5342
5343    // The union, in the order the spine runs: consecutive pieces meet on the
5344    // *same* placed profile, which is the coincident-face case the boolean
5345    // resolves by identifying it rather than by probing across it.
5346    let mut history = History::new();
5347    let mut shape = pieces[0].clone();
5348    for piece in &pieces[1..] {
5349        shape = ogeom_bool::fuse(model, &shape, piece, tol)?.shape;
5350    }
5351    history.generate(spine, shape.clone());
5352    history.generate(profile, shape.clone());
5353    Ok(Built::new(shape, history))
5354}
5355
5356/// The profile as a face.
5357///
5358/// A closed profile bounds its own area. An open one does not, and there is
5359/// exactly one honest way to close it: against the plane the spine was given
5360/// as a face *in*, which is what a face spine says to do and a wire spine has
5361/// no answer for. The closing segment runs between the profile's two ends, and
5362/// both have to be on that plane or the profile does not reach it.
5363fn profile_face(
5364    model: &mut Model,
5365    profile: &Shape,
5366    profile_normal: Vector,
5367    against: Option<Plane>,
5368    tol: Tolerances,
5369) -> OgeomResult<Shape> {
5370    if model.kind_of(profile)? == ShapeType::Face {
5371        return Ok(profile.clone());
5372    }
5373    if model.kind_of(profile)? != ShapeType::Wire {
5374        ogeom_bail!(Construction, "a profile is a wire or a face");
5375    }
5376    let mut edges = ogeom_topo::explore(model, profile, Filter::OfType(ShapeType::Edge))?;
5377    let closed = ogeom_algo::is_wire_closed(model, profile, tol)?;
5378    if !closed {
5379        let Some(plane) = against else {
5380            ogeom_bail!(
5381                Construction,
5382                "an open profile sweeps a shell, not a volume; give the spine \
5383                 as a planar face for its plane to close the profile against, \
5384                 or close the profile itself"
5385            );
5386        };
5387        let [(from, v0), (to, v1)] = wire_ends(model, profile, tol)?;
5388        for end in [from, to] {
5389            if plane.signed_distance_to(end).abs() > tol.confusion() * 1e2 {
5390                ogeom_bail!(
5391                    Construction,
5392                    "an open profile is closed against the spine face's own \
5393                     plane, and this one does not reach it"
5394                );
5395            }
5396        }
5397        // Built on the profile's *own* end vertices, so the closed ring is a
5398        // wire rather than edges that merely touch.
5399        let line = LineCurve::new(ogeom_math::Axis {
5400            location: from,
5401            direction: Direction::new(to - from, tol)?,
5402        });
5403        edges.push(
5404            make_edge_between(model, line.into(), (0.0, from.distance(to)), &v0, &v1, tol)?.shape,
5405        );
5406    }
5407    let ordered = ogeom_algo::order_edges(model, &edges, tol)?;
5408    let mut bound = ogeom_math::Aabb::EMPTY;
5409    for edge in &ordered {
5410        bound = bound.union(&ogeom_algo::shape_bounds(model, edge, tol)?);
5411    }
5412    let Some(centre) = bound.centre() else {
5413        ogeom_bail!(Construction, "a profile with no extent sweeps nothing");
5414    };
5415    let reach = bound.diagonal().mul_add(2.0, 1.0);
5416    let plane = Plane::through(centre, Direction::new(profile_normal, tol)?);
5417    let surface = PlaneSurface::over(plane, (-reach, reach), (-reach, reach))?;
5418    Ok(make_face_with_pcurves(model, surface.into(), &[ordered], tol)?.shape)
5419}
5420
5421/// Where an open wire begins and ends: the point, and the vertex there.
5422fn wire_ends(model: &Model, wire: &Shape, tol: Tolerances) -> OgeomResult<[(Point, Shape); 2]> {
5423    let mut counts: Vec<(Point, Shape, usize)> = Vec::new();
5424    for edge in explore(model, wire, Filter::OfType(ShapeType::Edge))? {
5425        for v in explore(model, &edge, Filter::OfType(ShapeType::Vertex))? {
5426            let Some(data) = model.node(&v).and_then(|n| n.data().as_vertex()) else {
5427                continue;
5428            };
5429            let at = v.transform(model.datums())?.apply(data.point);
5430            match counts
5431                .iter_mut()
5432                .find(|(p, _, _)| p.distance(at) <= tol.confusion() * 10.0)
5433            {
5434                Some((_, _, n)) => *n += 1,
5435                None => counts.push((at, v.clone(), 1)),
5436            }
5437        }
5438    }
5439    let free: Vec<(Point, Shape)> = counts
5440        .into_iter()
5441        .filter(|(_, _, n)| *n == 1)
5442        .map(|(p, v, _)| (p, v))
5443        .collect();
5444    if free.len() != 2 {
5445        ogeom_bail!(
5446            Construction,
5447            "an open profile has exactly two ends; this one has {}",
5448            free.len()
5449        );
5450    }
5451    let mut ends = free.into_iter();
5452    let (Some(a), Some(b)) = (ends.next(), ends.next()) else {
5453        ogeom_bail!(Construction, "the profile lost an end between checks");
5454    };
5455    Ok([a, b])
5456}
5457
5458/// The spine, edge by edge, in the order the wire runs it.
5459fn spine_stations(model: &Model, wire: &Shape, tol: Tolerances) -> OgeomResult<Vec<Station>> {
5460    let mut out = Vec::new();
5461    for edge in explore(model, wire, Filter::OfType(ShapeType::Edge))? {
5462        let Some(data) = model.node(&edge).and_then(|n| n.data().as_edge()) else {
5463            ogeom_bail!(Construction, "a spine edge is not in this model");
5464        };
5465        let Some(EdgeRepr::Curve3d { curve, range, .. }) = data.curve3d() else {
5466            ogeom_bail!(Construction, "a spine edge with no curve runs nowhere");
5467        };
5468        let Some(geometry) = model.geometry().curve(*curve) else {
5469            ogeom_bail!(Dangling, "curve is not in this model");
5470        };
5471        let placed = geometry.transformed(&edge.transform(model.datums())?, tol)?;
5472        let reversed = edge.orientation() == ogeom_topo::Orientation::Reversed;
5473        let (t0, t1) = if reversed {
5474            (range.1, range.0)
5475        } else {
5476            (range.0, range.1)
5477        };
5478        let sign = if reversed { -1.0 } else { 1.0 };
5479        let unit = |t: f64| -> OgeomResult<Vector> {
5480            let d = placed.d1_at(t, tol)? * sign;
5481            if d.magnitude() <= tol.confusion() {
5482                ogeom_bail!(Construction, "a spine edge has no direction at {t}");
5483            }
5484            Ok(d / d.magnitude())
5485        };
5486        let station = match &placed {
5487            Curve::Line(_) => Station {
5488                from: placed.point_at(t0, tol)?,
5489                to: placed.point_at(t1, tol)?,
5490                tangent_in: unit(t0)?,
5491                tangent_out: unit(t1)?,
5492                turn: None,
5493            },
5494            Curve::Circle(c) => {
5495                let circle = c.circle();
5496                let swept = (range.1 - range.0).abs();
5497                let axis = ogeom_math::Axis {
5498                    location: circle.centre(),
5499                    direction: if reversed {
5500                        -circle.frame().z()
5501                    } else {
5502                        circle.frame().z()
5503                    },
5504                };
5505                Station {
5506                    from: placed.point_at(t0, tol)?,
5507                    to: placed.point_at(t1, tol)?,
5508                    tangent_in: unit(t0)?,
5509                    tangent_out: unit(t1)?,
5510                    turn: Some((axis, swept)),
5511                }
5512            }
5513            other => ogeom_bail!(
5514                Construction,
5515                "a spine runs on straight and circular edges; a {:?} sweeps a \
5516                 surface this construction does not have",
5517                other.kind()
5518            ),
5519        };
5520        out.push(station);
5521    }
5522    Ok(out)
5523}
5524
5525/// The spine's own normal, and the check that it has one.
5526///
5527/// Taken from the first turn the spine makes (a corner or an arc) because
5528/// that is exact, and then measured against every station: a spine that
5529/// leaves its own plane has no square profile to carry, and says so here
5530/// rather than by producing a shape nobody asked for.
5531fn spine_normal(stations: &[Station], tol: Tolerances) -> OgeomResult<Direction> {
5532    let mut best: Option<(f64, Vector)> = None;
5533    let mut consider = |a: Vector, b: Vector| {
5534        let cross = a.cross(b);
5535        let magnitude = cross.magnitude();
5536        if magnitude > best.map_or(tol.angular(), |(m, _)| m) {
5537            best = Some((magnitude, cross / magnitude));
5538        }
5539    };
5540    for (index, station) in stations.iter().enumerate() {
5541        consider(station.tangent_in, station.tangent_out);
5542        if index + 1 < stations.len() {
5543            consider(station.tangent_out, stations[index + 1].tangent_in);
5544        }
5545    }
5546    if stations.len() > 1 {
5547        consider(
5548            stations[stations.len() - 1].tangent_out,
5549            stations[0].tangent_in,
5550        );
5551    }
5552    let Some((_, normal)) = best else {
5553        ogeom_bail!(
5554            Construction,
5555            "a spine that never turns has no plane of its own; give the \
5556             profile's own orientation a spine with at least one corner or arc"
5557        );
5558    };
5559    for station in stations {
5560        for tangent in [station.tangent_in, station.tangent_out] {
5561            if tangent.dot(normal).abs() > tol.angular() {
5562                ogeom_bail!(
5563                    Construction,
5564                    "the spine leaves its own plane; an evolved sweep runs a \
5565                     planar spine"
5566                );
5567            }
5568        }
5569        if let Some((axis, _)) = &station.turn
5570            && axis.direction.vector().cross(normal).magnitude() > tol.angular()
5571        {
5572            ogeom_bail!(
5573                Construction,
5574                "a spine arc turns about an axis off the spine's own normal"
5575            );
5576        }
5577    }
5578    Direction::new(normal, tol)
5579}
5580
5581/// The profile's plane: a point on it and its normal.
5582fn profile_plane(model: &Model, profile: &Shape, tol: Tolerances) -> OgeomResult<(Point, Vector)> {
5583    let mut points: Vec<Point> = Vec::new();
5584    for edge in explore(model, profile, Filter::OfType(ShapeType::Edge))? {
5585        let Some(data) = model.node(&edge).and_then(|n| n.data().as_edge()) else {
5586            continue;
5587        };
5588        let Some(EdgeRepr::Curve3d { curve, range, .. }) = data.curve3d() else {
5589            continue;
5590        };
5591        let Some(geometry) = model.geometry().curve(*curve) else {
5592            ogeom_bail!(Dangling, "curve is not in this model");
5593        };
5594        let placed = geometry.transformed(&edge.transform(model.datums())?, tol)?;
5595        for i in 0..=8 {
5596            let t = range.0 + (range.1 - range.0) * f64::from(i) / 8.0;
5597            points.push(placed.point_at(t, tol)?);
5598        }
5599    }
5600    if points.len() < 3 {
5601        ogeom_bail!(Construction, "a profile needs an extent to sweep");
5602    }
5603    let origin = points[0];
5604    // The widest cross product among the sampled offsets: the plane's normal,
5605    // taken where it is best conditioned rather than from the first three
5606    // points that happen to be there.
5607    let mut best: Option<(f64, Vector)> = None;
5608    for (i, a) in points.iter().enumerate() {
5609        for b in points.iter().skip(i + 1) {
5610            let cross = (*a - origin).cross(*b - origin);
5611            let magnitude = cross.magnitude();
5612            if magnitude > best.map_or(tol.confusion(), |(m, _)| m) {
5613                best = Some((magnitude, cross / magnitude));
5614            }
5615        }
5616    }
5617    let Some((_, normal)) = best else {
5618        ogeom_bail!(Construction, "a profile with no area has no plane");
5619    };
5620    for p in &points {
5621        if (*p - origin).dot(normal).abs() > tol.confusion() * 1e2 {
5622            ogeom_bail!(Construction, "the profile is not planar");
5623        }
5624    }
5625    Ok((origin, normal))
5626}
5627
5628/// The rigid motion that carries the profile from the spine's start to a
5629/// station: a turn about the spine's normal, then a translation.
5630fn station_transform(
5631    reference: (Point, Vector),
5632    at: Point,
5633    tangent: Vector,
5634    normal: Direction,
5635    tol: Tolerances,
5636) -> OgeomResult<Transform> {
5637    let (origin, from) = reference;
5638    let n = normal.vector();
5639    let angle = from.cross(tangent).dot(n).atan2(from.dot(tangent));
5640    let turn = if angle.abs() <= tol.angular() {
5641        Transform::IDENTITY
5642    } else {
5643        Transform::rotation(
5644            ogeom_math::Axis {
5645                location: origin,
5646                direction: normal,
5647            },
5648            angle,
5649        )
5650    };
5651    Ok(Transform::translation(at - origin) * turn)
5652}
5653
5654/// The wedge a corner adds: the profile turned about the corner, through
5655/// exactly the angle the spine turns there.
5656///
5657/// `None` where the spine does not turn; two edges meeting smoothly leave no
5658/// wedge to fill.
5659#[allow(clippy::too_many_arguments)]
5660fn corner_piece(
5661    model: &mut Model,
5662    profile: &Shape,
5663    reference: (Point, Vector),
5664    corner: Point,
5665    incoming: Vector,
5666    outgoing: Vector,
5667    normal: Direction,
5668    tol: Tolerances,
5669) -> OgeomResult<Option<Shape>> {
5670    let n = normal.vector();
5671    let angle = incoming
5672        .cross(outgoing)
5673        .dot(n)
5674        .atan2(incoming.dot(outgoing));
5675    if angle.abs() <= tol.angular() {
5676        return Ok(None);
5677    }
5678    let placed = ogeom_algo::transformed(
5679        model,
5680        profile,
5681        station_transform(reference, corner, incoming, normal, tol)?,
5682    )?
5683    .shape;
5684    let axis = ogeom_math::Axis {
5685        location: corner,
5686        direction: if angle > 0.0 { normal } else { -normal },
5687    };
5688    let turned = ogeom_algo::make_revolution(model, &placed, axis, angle.abs(), tol);
5689    match turned {
5690        Ok(built) => Ok(Some(built.shape)),
5691        Err(_) => ogeom_bail!(
5692            NotDone,
5693            "the profile straddles the spine at a corner, so turning it about \
5694             that corner sweeps it through itself; there is no revolution for \
5695             that wedge"
5696        ),
5697    }
5698}