1use ogeom_algo::Built;
15use ogeom_core::{OgeomResult, Tolerances, ogeom_bail};
16use ogeom_geom::{PlaneSurface, Surface as _, SurfaceGeometry};
17use ogeom_math::{Direction, Frame, Plane, Point, Transform, Vector};
18use ogeom_topo::{Model, NodeData, Shape, ShapeType, TShapeId};
19
20use crate::shape::rebuilt;
21
22pub fn apply_draft(
43 model: &mut Model,
44 solid: &Shape,
45 faces: &[Shape],
46 neutral: Plane,
47 pull: Direction,
48 angle: f64,
49 tol: Tolerances,
50) -> OgeomResult<Built> {
51 if !angle.is_finite() || angle.abs() >= core::f64::consts::FRAC_PI_2 {
52 ogeom_bail!(
53 Construction,
54 "a draft of {angle} radians turns the face past its own plane"
55 );
56 }
57 let (canonical, mapped, prefix) = crate::shape::canonical_input(model, solid, faces, tol)?;
58 if let Some(prefix) = prefix {
59 let mut out = apply_draft(model, &canonical, &mapped, neutral, pull, angle, tol)?;
60 out.history = prefix.then(&out.history);
61 return Ok(out);
62 }
63 if faces.is_empty() {
64 ogeom_bail!(Construction, "a draft of no faces drafts nothing");
65 }
66 let own: Vec<Shape> = {
70 let mut seen: Vec<Shape> = Vec::new();
71 for f in ogeom_topo::explore(model, solid, ogeom_topo::Filter::OfType(ShapeType::Face))? {
72 if !seen.iter().any(|s| s.node() == f.node()) {
73 seen.push(f);
74 }
75 }
76 seen
77 };
78
79 let mut turned: Vec<(TShapeId, SurfaceGeometry)> = Vec::with_capacity(faces.len());
83 for face in faces {
84 let Some(used) = own.iter().find(|f| f.node() == face.node()).cloned() else {
85 ogeom_bail!(Construction, "a drafted face is not a face of the solid");
86 };
87 let face = &used;
88 let Some(NodeData::Face(data)) = model.node(face).map(|n| n.data().clone()) else {
89 ogeom_bail!(Construction, "expected a face");
90 };
91 let Some(surface) = model.geometry().surface(data.surface) else {
92 ogeom_bail!(Dangling, "face refers to a surface not in this model");
93 };
94 let sign = outward_sign(model, solid, face, surface, tol)?;
99 let sign_of = |_: &Shape| sign;
100 let axial = |frame: Frame| -> bool {
103 (frame.z().vector().dot(neutral.normal().vector()).abs() - 1.0).abs()
104 <= tol.angular().max(1e-9)
105 };
106 match surface {
107 SurfaceGeometry::Cylinder(c) if axial(c.cylinder().frame()) => {
108 let cylinder = c.cylinder();
109 let (_, (v0, v1)) = surface.domain();
110 turned.push((
111 face.node(),
112 revolved_draft(
113 cylinder.frame(),
114 cylinder.radius(),
115 0.0,
116 (v0, v1),
117 sign_of(face),
118 neutral,
119 pull,
120 angle,
121 tol,
122 )?,
123 ));
124 continue;
125 }
126 SurfaceGeometry::Cone(co) if axial(co.cone().frame()) => {
127 let cone = co.cone();
128 let (_, (v0, v1)) = surface.domain();
129 turned.push((
130 face.node(),
131 revolved_draft(
132 cone.frame(),
133 cone.reference_radius(),
134 cone.half_angle(),
135 (v0, v1),
136 sign_of(face),
137 neutral,
138 pull,
139 angle,
140 tol,
141 )?,
142 ));
143 continue;
144 }
145 SurfaceGeometry::Extrusion(e) => {
146 turned.push((
147 face.node(),
148 extruded_draft(
149 e,
150 surface.domain(),
151 sign_of(face),
152 neutral,
153 pull,
154 angle,
155 tol,
156 )?,
157 ));
158 continue;
159 }
160 SurfaceGeometry::Plane(_) => {}
161 _ => {
166 turned.push((
167 face.node(),
168 general_draft(model, face, surface, sign, neutral, pull, angle, tol)?,
169 ));
170 continue;
171 }
172 }
173 let SurfaceGeometry::Plane(p) = surface else {
174 unreachable!("the match above let only planes through");
175 };
176 let plane = p.plane();
177 let ((u0, u1), (v0, v1)) = surface.domain();
178 let outward = plane.normal().vector() * sign;
181
182 let along = plane.normal().vector().cross(neutral.normal().vector());
184 let magnitude = along.magnitude();
185 if magnitude <= tol.angular() {
186 ogeom_bail!(
187 Construction,
188 "a face parallel to the neutral plane has no line to turn \
189 about"
190 );
191 }
192 let along = along / magnitude;
193 let hinge = meet(plane, neutral, along, tol)?;
194
195 let axis = ogeom_math::Axis::new(hinge, Direction::new(along, tol)?);
200 let mut candidates = Vec::with_capacity(2);
201 for sense in [1.0, -1.0] {
202 let turn = Transform::rotation(axis, angle.abs() * sense);
203 candidates.push((sense, turn.apply_vector(outward).dot(pull.vector())));
204 }
205 let leaning = candidates
206 .iter()
207 .copied()
208 .max_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(core::cmp::Ordering::Equal))
209 .map_or(1.0, |(sense, _)| sense);
210 let turn = Transform::rotation(axis, angle * leaning);
211 let moved_normal = Direction::new(turn.apply_vector(plane.normal().vector()), tol)?;
212 let tilted = Plane::new(Frame::new(
213 hinge,
214 moved_normal,
215 Direction::new(along, tol)?,
216 tol,
217 )?);
218 let ((u0, u1), (v0, v1)) = if [u0, u1, v0, v1].iter().all(|x| x.is_finite()) {
222 let (mut lo, mut hi) = (
223 (f64::INFINITY, f64::INFINITY),
224 (f64::NEG_INFINITY, f64::NEG_INFINITY),
225 );
226 for (u, v) in [(u0, v0), (u1, v0), (u1, v1), (u0, v1)] {
227 let corner = tilted
228 .frame()
229 .to_local(turn.apply(surface.point_at(u, v, tol)?));
230 lo = (lo.0.min(corner.x), lo.1.min(corner.y));
231 hi = (hi.0.max(corner.x), hi.1.max(corner.y));
232 }
233 ((lo.0, hi.0), (lo.1, hi.1))
234 } else {
235 ((u0, u1), (v0, v1))
236 };
237 let grow = (u1 - u0).abs().max((v1 - v0).abs()).mul_add(0.5, 1.0) * angle.abs().tan()
240 + tol.confusion();
241 turned.push((
242 face.node(),
243 PlaneSurface::over(tilted, (u0 - grow, u1 + grow), (v0 - grow, v1 + grow))?.into(),
244 ));
245 }
246
247 rebuilt(
248 model,
249 solid,
250 &|_| 0.0,
251 &|face| {
252 turned
253 .iter()
254 .find(|(node, _)| *node == face.node())
255 .map(|(_, surface)| surface.clone())
256 },
257 tol,
258 )
259}
260
261#[allow(clippy::too_many_arguments, reason = "one construction, all its data")]
265fn revolved_draft(
266 frame: Frame,
267 reference_radius: f64,
268 half_angle: f64,
269 window: (f64, f64),
270 sign: f64,
271 neutral: Plane,
272 pull: Direction,
273 angle: f64,
274 tol: Tolerances,
275) -> OgeomResult<SurfaceGeometry> {
276 use ogeom_geom::ConeSurface;
277
278 let axis_dir = frame.z().vector();
279 let along = axis_dir.dot(neutral.normal().vector());
280 if (along.abs() - 1.0).abs() > tol.angular().max(1e-9) {
281 ogeom_bail!(
282 Construction,
283 "a wall of revolution drafts about a neutral plane square to \
284 its axis; the oblique neutral needs the general machinery (see \
285 docs/CAPABILITIES.md, offset.draft)"
286 );
287 }
288 let height = -neutral.signed_distance_to(frame.origin()) * along.signum();
291 let neutral_point = frame.origin() + axis_dir * height;
292 let neutral_radius = half_angle.tan().mul_add(height, reference_radius);
293 if neutral_radius <= tol.confusion() {
294 ogeom_bail!(
295 Construction,
296 "the wall has no radius left at the neutral plane to hold"
297 );
298 }
299 let hinge_frame = Frame::new(neutral_point, frame.z(), frame.x(), tol)?;
300
301 let mut best: Option<(f64, f64)> = None;
305 for sense in [1.0_f64, -1.0] {
306 let probe = half_angle + angle.abs() * sense;
309 let candidate = half_angle + angle * sense;
310 if probe.abs() <= tol.angular()
311 || probe.abs() >= core::f64::consts::FRAC_PI_2 - tol.angular()
312 || candidate.abs() <= tol.angular()
313 || candidate.abs() >= core::f64::consts::FRAC_PI_2 - tol.angular()
314 {
315 continue;
316 }
317 let cone = ogeom_math::Cone::new(hinge_frame, neutral_radius, probe, tol)?;
318 let surface: SurfaceGeometry = ConeSurface::new(cone, (-1.0, 1.0))?.into();
319 let (du, dv) = surface.d1_at(0.0, 1.0, tol)?;
320 let n = du.cross(dv);
321 let outward = n / n.magnitude() * sign;
322 let lean = outward.dot(pull.vector());
323 if best.as_ref().is_none_or(|(_, held)| lean > *held) {
324 best = Some((candidate, lean));
325 }
326 }
327 let Some((leaned, _)) = best else {
328 ogeom_bail!(
329 Construction,
330 "a draft of {angle} radians flattens the wall or swallows it"
331 );
332 };
333 let cone = ogeom_math::Cone::new(hinge_frame, neutral_radius, leaned, tol)?;
334
335 let shift = height;
338 let grow = (window.1 - window.0).abs().mul_add(0.1, 1.0);
339 let (w0, w1) = (window.0 - shift - grow, window.1 - shift + grow);
340 let apex_height = -neutral_radius / leaned.tan();
341 if apex_height > w0 && apex_height < w1 {
342 ogeom_bail!(
343 Construction,
344 "the draft swallows the drafted face's own apex"
345 );
346 }
347 Ok(ConeSurface::new(cone, (w0, w1))?.into())
348}
349
350#[allow(clippy::too_many_arguments, reason = "one construction, all its data")]
362fn extruded_draft(
363 extrusion: &ogeom_geom::ExtrusionSurface,
364 window: ((f64, f64), (f64, f64)),
365 sign: f64,
366 neutral: Plane,
367 pull: Direction,
368 angle: f64,
369 tol: Tolerances,
370) -> OgeomResult<SurfaceGeometry> {
371 use ogeom_geom::Curve3d as _;
372 let ((u0, u1), (v0, v1)) = window;
373 let d = extrusion.direction().vector();
374 let n = neutral.normal().vector();
375 let den = n.dot(d);
376 if den.abs() <= tol.angular() {
377 ogeom_bail!(
378 Construction,
379 "the neutral plane runs along the wall's rulings; there is no \
380 hinge to turn about"
381 );
382 }
383 let curve = extrusion.curve();
384 let o = neutral.origin().to_vector();
385 let height_at = |c: Point| n.dot(o - c.to_vector()) / den;
388 let hinge_tangent = |cd: Vector| cd - d * (n.dot(cd) / den);
389
390 let um = f64::midpoint(u0, u1);
395 let cm = curve.point_at(um, tol)?;
396 let cdm = curve.d1_at(um, tol)?;
397 let hinge_m = cm + d * height_at(cm);
398 let tangent_m = Direction::new(hinge_tangent(cdm), tol)?;
399 let outward = {
400 let nw = cdm.cross(d);
401 nw / nw.magnitude() * sign
402 };
403 let axis_m = ogeom_math::Axis::new(hinge_m, tangent_m);
404 let mut leaning = 1.0;
405 let mut best = f64::NEG_INFINITY;
406 for sense in [1.0_f64, -1.0] {
407 let turn = Transform::rotation(axis_m, angle.abs() * sense);
408 let lean = turn.apply_vector(outward).dot(pull.vector());
409 if lean > best {
410 best = lean;
411 leaning = sense;
412 }
413 }
414 let theta = angle * leaning;
415
416 let ruling_at = |u: f64| -> OgeomResult<(Point, Vector, f64)> {
419 let c = curve.point_at(u, tol)?;
420 let cd = curve.d1_at(u, tol)?;
421 let h = height_at(c);
422 let hinge = c + d * h;
423 let tangent = Direction::new(hinge_tangent(cd), tol)?;
424 let turn = Transform::rotation(ogeom_math::Axis::new(hinge, tangent), theta);
425 Ok((hinge, turn.apply_vector(d), h))
426 };
427 const ALONG: usize = 65;
428 #[allow(clippy::cast_precision_loss)]
429 let along_u = |x: f64| u0 + (u1 - u0) * x / ((ALONG - 1) as f64);
430 let mut base: Vec<(Point, Vector)> = Vec::with_capacity(ALONG);
431 let (mut s_lo, mut s_hi) = (f64::INFINITY, f64::NEG_INFINITY);
432 for i in 0..ALONG {
433 #[allow(clippy::cast_precision_loss)]
434 let (hinge, ruling, h) = ruling_at(along_u(i as f64))?;
435 base.push((hinge, ruling));
436 s_lo = s_lo.min(v0 - h);
437 s_hi = s_hi.max(v1 - h);
438 }
439 let grow = (u1 - u0).abs().max((v1 - v0).abs()).mul_add(0.5, 1.0) * angle.abs().tan()
444 + tol.confusion();
445 let (s_lo, s_hi) = (s_lo - grow, s_hi + grow);
446 struct Continuation {
450 hinge: Point,
451 d1: Vector,
452 d2: Vector,
453 ruling: Vector,
454 r1: Vector,
455 steps: usize,
456 }
457 let continuation = |front: bool| -> Continuation {
458 let (i0, i1, i2) = if front {
459 (0, 1, 2)
460 } else {
461 (ALONG - 1, ALONG - 2, ALONG - 3)
462 };
463 let d1 = base[i0].0 - base[i1].0;
464 let d2 = (base[i0].0 - base[i1].0) - (base[i1].0 - base[i2].0);
465 let step = d1.magnitude().max(tol.confusion());
470 let ruling = base[i0].1;
471 let lean = (ruling - d * ruling.dot(d)).magnitude() * s_lo.abs().max(s_hi.abs());
472 let steps = ((grow + lean) / step).ceil().max(2.0);
473 #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
474 let steps = (steps as usize).min(16);
475 Continuation {
476 hinge: base[i0].0,
477 d1,
478 d2,
479 ruling: base[i0].1,
480 r1: base[i0].1 - base[i1].1,
481 steps,
482 }
483 };
484 let (front, back) = (continuation(true), continuation(false));
485 let continued = |c: &Continuation, k: f64| -> (Point, Vector) {
486 (
487 c.hinge + c.d1 * k + c.d2 * (k * (k + 1.0) / 2.0),
488 c.ruling + c.r1 * k,
489 )
490 };
491 #[allow(clippy::cast_precision_loss)]
496 let (first, last) = (front.steps as f64, (front.steps + ALONG - 1) as f64);
497 let station = |x: f64| -> OgeomResult<(Point, Vector)> {
498 if x < first {
499 return Ok(continued(&front, first - x));
500 }
501 if x > last {
502 return Ok(continued(&back, x - last));
503 }
504 let (hinge, ruling, _) = ruling_at(along_u(x - first))?;
505 Ok((hinge, ruling))
506 };
507 let along_total = front.steps + ALONG + back.steps;
508 let stations: Vec<(Point, Vector)> = (0..along_total)
509 .map(|i| {
510 #[allow(clippy::cast_precision_loss)]
511 let x = i as f64;
512 match i.checked_sub(front.steps) {
513 Some(j) if j < ALONG => Ok(base[j]),
514 _ => station(x),
515 }
516 })
517 .collect::<OgeomResult<_>>()?;
518
519 for edge in [s_lo, s_hi] {
522 for pair in stations.windows(2) {
523 let ((h0, r0), (h1, r1)) = (pair[0], pair[1]);
524 let step = (h1 + r1 * edge) - (h0 + r0 * edge);
525 if step.dot(h1 - h0) <= 0.0 {
526 ogeom_bail!(
527 Construction,
528 "the draft folds the wall onto itself inside the drafted \
529 window; refused; see docs/CAPABILITIES.md, offset.draft"
530 );
531 }
532 }
533 }
534
535 const ACROSS: usize = 9;
540 #[allow(clippy::cast_precision_loss)]
545 let span = (along_total - 1) as f64;
546 let xs: Vec<f64> = (0..along_total)
547 .map(|i| {
548 #[allow(clippy::cast_precision_loss)]
549 let x = i as f64 / span;
550 x
551 })
552 .collect();
553 let ss: Vec<f64> = (0..ACROSS)
554 .map(|j| {
555 #[allow(clippy::cast_precision_loss)]
556 let s = j as f64 / ((ACROSS - 1) as f64);
557 s
558 })
559 .collect();
560 let fit_target = (tol.confusion() * 1e3).max(1e-4);
561 let fitted = ogeom_geom::fit::fit_surface_sampled(
562 |s, x| {
563 let x = x * span;
564 #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
565 let i = x.round() as usize;
566 #[allow(clippy::cast_precision_loss)]
567 let (hinge, ruling) = if (i as f64 - x).abs() <= 1e-12 * span && i < along_total {
568 stations[i]
569 } else {
570 station(x)?
571 };
572 Ok(hinge + ruling * (s_hi - (s_hi - s_lo) * s))
573 },
574 &ogeom_geom::fit::Sampling {
575 us: ss,
576 vs: xs,
577 between: (false, true),
578 closed_v: false,
579 most: 4096,
580 },
581 3,
582 fit_target,
583 tol,
584 )?;
585 if !fitted.met {
586 ogeom_bail!(
587 NotDone,
588 "the drafted wall's fit reached {} against a target of {fit_target}",
589 fitted.error
590 );
591 }
592 Ok(fitted.curve.into())
593}
594
595const HINGE_STATIONS: usize = 256;
597
598fn reaching(
605 surface: &SurfaceGeometry,
606 neutral: Plane,
607 tol: Tolerances,
608) -> OgeomResult<SurfaceGeometry> {
609 let ((ua, ub), (va, vb)) = surface.domain();
610 let apex = match surface {
611 SurfaceGeometry::Cylinder(_) => None,
612 SurfaceGeometry::Cone(c) => {
613 let cone = c.cone();
614 Some(-cone.reference_radius() / cone.half_angle().tan())
615 }
616 _ => return Ok(surface.clone()),
617 };
618 let n = neutral.normal().vector();
619 let o = neutral.origin().to_vector();
620 let (mut lo, mut hi) = (va, vb);
621 for k in 0..HINGE_STATIONS {
622 #[allow(clippy::cast_precision_loss)]
623 let u = (ub - ua).mul_add(k as f64 / HINGE_STATIONS as f64, ua);
624 let foot = surface.point_at(u, 0.0, tol)?;
625 let ruling = surface.point_at(u, 1.0, tol)? - foot;
626 let along = n.dot(ruling);
627 if along.abs() <= tol.angular() {
628 continue;
629 }
630 let v = n.dot(o - foot.to_vector()) / along;
631 (lo, hi) = (lo.min(v), hi.max(v));
632 }
633 let margin = (hi - lo) * 0.1;
634 let (lo, hi) = (lo - margin, hi + margin);
635 if let Some(apex) = apex
637 && ((va >= apex && lo < apex) || (vb <= apex && hi > apex))
638 {
639 return Ok(surface.clone());
640 }
641 Ok(match surface {
642 SurfaceGeometry::Cylinder(c) => {
643 ogeom_geom::CylinderSurface::new(c.cylinder(), (lo, hi))?.into()
644 }
645 SurfaceGeometry::Cone(c) => ogeom_geom::ConeSurface::new(c.cone(), (lo, hi))?.into(),
646 _ => surface.clone(),
647 })
648}
649
650fn crossing_segments(
657 mesh: &ogeom_topo::Triangulation,
658 neutral: Plane,
659 on: f64,
660) -> Vec<[((f64, f64), Point); 2]> {
661 let side: Vec<f64> = mesh
662 .positions
663 .iter()
664 .map(|p| neutral.signed_distance_to(*p))
665 .collect();
666 let mut segments: Vec<[((f64, f64), Point); 2]> = Vec::new();
667 for t in &mesh.triangles {
668 let mut ends: Vec<((f64, f64), Point)> = Vec::with_capacity(2);
669 for &corner in t {
670 let i = corner as usize;
671 if side[i].abs() <= on {
672 ends.push((mesh.parameters[i], mesh.positions[i]));
673 }
674 }
675 for k in 0..3 {
676 let (i, j) = (t[k] as usize, t[(k + 1) % 3] as usize);
677 let (a, b) = (side[i], side[j]);
678 if a.abs() <= on || b.abs() <= on || (a < 0.0) == (b < 0.0) {
679 continue;
680 }
681 let f = a / (a - b);
682 let (pa, pb) = (mesh.parameters[i], mesh.parameters[j]);
683 let (qa, qb) = (mesh.positions[i], mesh.positions[j]);
684 ends.push((
685 (pa.0 + (pb.0 - pa.0) * f, pa.1 + (pb.1 - pa.1) * f),
686 qa + (qb - qa) * f,
687 ));
688 }
689 ends.dedup_by(|a, b| a.1.distance(b.1) <= on);
690 if ends.len() == 2 && ends[0].1.distance(ends[1].1) > on {
691 segments.push([ends[0], ends[1]]);
692 }
693 }
694 segments
695}
696
697#[allow(clippy::too_many_arguments, reason = "one construction, all its data")]
706fn general_draft(
707 model: &Model,
708 face: &Shape,
709 surface: &SurfaceGeometry,
710 sign: f64,
711 neutral: Plane,
712 pull: Direction,
713 angle: f64,
714 tol: Tolerances,
715) -> OgeomResult<SurfaceGeometry> {
716 let n = neutral.normal().vector();
717 let mesh = ogeom_mesh::triangulate_face(
718 model,
719 face,
720 ogeom_mesh::Deflection {
721 chord: 0.05,
722 angular: 0.2,
723 ..ogeom_mesh::Deflection::default()
724 },
725 tol,
726 )?;
727 let extent = {
728 let b = mesh
729 .positions
730 .iter()
731 .fold(ogeom_math::Aabb::EMPTY, |acc, p| acc.with_point(*p));
732 match (b.low(), b.high()) {
733 (Some(lo), Some(hi)) => (hi - lo).magnitude(),
734 _ => ogeom_bail!(Construction, "the drafted face has no extent"),
735 }
736 };
737 if extent <= tol.confusion() {
738 ogeom_bail!(Construction, "the drafted face has no extent");
739 }
740
741 let on = tol.confusion() * 10.0;
747 let mut segments = crossing_segments(&mesh, neutral, on);
748 let beyond = segments.is_empty();
749 let reached;
752 let surface: &SurfaceGeometry = if beyond {
753 reached = reaching(surface, neutral, tol)?;
754 &reached
755 } else {
756 surface
757 };
758 if beyond {
759 let mut scratch = Model::new();
760 let whole = ogeom_algo::make_natural_face(&mut scratch, surface.clone())?.shape;
761 let whole_mesh = ogeom_mesh::triangulate_face(
762 &scratch,
763 &whole,
764 ogeom_mesh::Deflection {
765 chord: 0.05,
766 angular: 0.2,
767 ..ogeom_mesh::Deflection::default()
768 },
769 tol,
770 )?;
771 segments = crossing_segments(&whole_mesh, neutral, on);
772 }
773 if segments.is_empty() {
774 ogeom_bail!(
775 Construction,
776 "the neutral plane meets neither the drafted face nor its surface; \
777 there is no hinge to turn about"
778 );
779 }
780 let same = |a: Point, b: Point| a.distance(b) <= on;
786 let mut chain: Vec<((f64, f64), Point)> = vec![segments[0][0], segments[0][1]];
787 let mut used = vec![false; segments.len()];
788 used[0] = true;
789 loop {
790 let mut grew = false;
791 for (k, seg) in segments.iter().enumerate() {
792 if used[k] {
793 continue;
794 }
795 let tail = chain[chain.len() - 1].1;
799 let head = chain[0].1;
800 if (same(seg[0].1, tail) && same(seg[1].1, head))
803 || (same(seg[1].1, tail) && same(seg[0].1, head))
804 {
805 chain.push(chain[0]);
806 used[k] = true;
807 grew = true;
808 break;
809 }
810 let covered = |p: Point| chain.iter().any(|c| same(c.1, p));
811 if covered(seg[0].1) && covered(seg[1].1) {
812 used[k] = true;
813 continue;
814 }
815 if same(seg[0].1, tail) {
816 chain.push(seg[1]);
817 } else if same(seg[1].1, tail) {
818 chain.push(seg[0]);
819 } else if same(seg[0].1, head) {
820 chain.insert(0, seg[1]);
821 } else if same(seg[1].1, head) {
822 chain.insert(0, seg[0]);
823 } else {
824 continue;
825 }
826 used[k] = true;
827 grew = true;
828 }
829 if !grew {
830 break;
831 }
832 }
833 if used.iter().any(|u| !u) {
834 let what = if beyond {
839 "the drafted face's surface"
840 } else {
841 "the drafted face"
842 };
843 ogeom_bail!(
844 Construction,
845 "the neutral plane crosses {what} more than once; there is no one \
846 hinge to turn about"
847 );
848 }
849 let closed = chain.len() > 3 && same(chain[0].1, chain[chain.len() - 1].1);
850 if closed {
851 chain.pop();
852 }
853 if chain.len() < 2 {
854 ogeom_bail!(
855 Construction,
856 "the neutral plane touches the drafted face at a point; there is \
857 no hinge to turn about"
858 );
859 }
860 let chain: Vec<(f64, f64)> = chain.into_iter().map(|c| c.0).collect();
861
862 let ((ua, ub), (va, vb)) = surface.domain();
868 let period = (
872 (surface.is_periodic_u() || surface.is_closed_u(tol)).then_some(ub - ua),
873 (surface.is_periodic_v() || surface.is_closed_v(tol)).then_some(vb - va),
874 );
875 let short = |a: f64, b: f64, period: Option<f64>| -> f64 {
876 let d = b - a;
877 match period {
878 Some(p) if d.abs() > p * 0.5 => d - p * d.signum(),
879 _ => d,
880 }
881 };
882 let chain: Vec<(f64, f64)> = {
883 let chain: Vec<(f64, f64)> = if closed {
887 let n = chain.len();
891 let mut exact: Option<(usize, (f64, f64))> = None;
892 for i in 0..n {
893 let (a, b) = (chain[i], chain[(i + 1) % n]);
894 let da = short(ua, a.0, period.0);
898 let db = da + short(a.0, b.0, period.0);
899 if da == 0.0 {
900 exact = Some((i, a));
901 break;
902 }
903 if (da < 0.0) != (db < 0.0) && (da - db).abs() > 0.0 {
904 let f = da / (da - db);
905 let dv = short(a.1, b.1, period.1);
906 exact = Some((i + 1, (ua, a.1 + dv * f)));
907 break;
908 }
909 }
910 let (start, inserted) = exact.unwrap_or_else(|| {
911 let mut best = (0usize, f64::INFINITY);
912 for (i, c) in chain.iter().enumerate() {
913 let d = short(ua, c.0, period.0).abs();
914 if d < best.1 {
915 best = (i, d);
916 }
917 }
918 (best.0, chain[best.0])
919 });
920 let mut rotated: Vec<(f64, f64)> = Vec::with_capacity(n + 1);
921 rotated.push(inserted);
922 for k in 0..n {
923 let c = chain[(start + k) % n];
924 if rotated.len() == 1 && c == inserted {
925 continue;
926 }
927 rotated.push(c);
928 }
929 rotated
930 } else {
931 chain
932 };
933 let pairs = if closed { chain.len() } else { chain.len() - 1 };
937 let mut lengths = Vec::with_capacity(pairs);
938 let mut total = 0.0;
939 for i in 0..pairs {
940 let (a, b) = (chain[i], chain[(i + 1) % chain.len()]);
941 let step = surface
942 .point_at(a.0, a.1, tol)?
943 .distance(surface.point_at(b.0, b.1, tol)?);
944 lengths.push(step);
945 total += step;
946 }
947 let count = if closed {
948 HINGE_STATIONS
949 } else {
950 HINGE_STATIONS + 1
951 };
952 let mut dense = Vec::with_capacity(count);
953 let (mut pair, mut walked) = (0usize, 0.0_f64);
954 for k in 0..count {
955 #[allow(clippy::cast_precision_loss)]
956 let target = total * k as f64 / HINGE_STATIONS as f64;
957 while pair + 1 < pairs && walked + lengths[pair] < target {
958 walked += lengths[pair];
959 pair += 1;
960 }
961 let (a, b) = (chain[pair], chain[(pair + 1) % chain.len()]);
962 let (du, dv) = (short(a.0, b.0, period.0), short(a.1, b.1, period.1));
963 let f = if lengths[pair] > 0.0 {
964 ((target - walked) / lengths[pair]).clamp(0.0, 1.0)
965 } else {
966 0.0
967 };
968 dense.push((a.0 + du * f, a.1 + dv * f));
969 }
970 dense
971 };
972
973 let on_hinge = |(mut u, mut v): (f64, f64),
980 pinned: bool|
981 -> OgeomResult<((f64, f64), Point, Vector, Vector)> {
982 for _ in 0..8 {
983 let p = surface.point_at(u, v, tol)?;
984 let f = neutral.signed_distance_to(p);
985 if f.abs() <= tol.confusion() * 1e-2 {
986 break;
987 }
988 let (du, dv) = surface.d1_at(u, v, tol)?;
989 let g = (if pinned { 0.0 } else { n.dot(du) }, n.dot(dv));
990 let g2 = g.0 * g.0 + g.1 * g.1;
991 if g2 <= 0.0 {
992 break;
993 }
994 u -= f * g.0 / g2;
995 v -= f * g.1 / g2;
996 }
997 let p = surface.point_at(u, v, tol)?;
998 let (du, dv) = surface.d1_at(u, v, tol)?;
999 let raw = du.cross(dv);
1000 if raw.magnitude() <= tol.confusion() {
1001 ogeom_bail!(Construction, "the drafted face has no normal on its hinge");
1002 }
1003 let outward = raw / raw.magnitude() * sign;
1004 let t = outward.cross(n);
1006 if t.magnitude() <= tol.angular() {
1007 ogeom_bail!(
1008 Construction,
1009 "the neutral plane is tangent to the drafted face; there is no \
1010 hinge to turn about"
1011 );
1012 }
1013 Ok(((u, v), p, t / t.magnitude(), outward))
1014 };
1015 let mut feet: Vec<(f64, f64)> = Vec::with_capacity(chain.len());
1016 let mut hinges: Vec<Point> = Vec::with_capacity(chain.len());
1017 let mut tangents: Vec<Vector> = Vec::with_capacity(chain.len());
1018 let mut outwards: Vec<Vector> = Vec::with_capacity(chain.len());
1019 for (k, &at) in chain.iter().enumerate() {
1020 let (foot, p, mut t, outward) = on_hinge(at, closed && k == 0)?;
1023 let next = chain[(k + 1) % chain.len()];
1025 let prev = chain[(k + chain.len() - 1) % chain.len()];
1026 let ahead =
1027 surface.point_at(next.0, next.1, tol)? - surface.point_at(prev.0, prev.1, tol)?;
1028 if t.dot(ahead) < 0.0 {
1029 t = -t;
1030 }
1031 feet.push(foot);
1032 hinges.push(p);
1033 tangents.push(t);
1034 outwards.push(outward);
1035 }
1036
1037 let m = hinges.len() / 2;
1041 let axis_m = ogeom_math::Axis::new(hinges[m], Direction::new(tangents[m], tol)?);
1042 let mut leaning = 1.0;
1043 let mut best = f64::NEG_INFINITY;
1044 for sense in [1.0_f64, -1.0] {
1045 let turn = Transform::rotation(axis_m, angle.abs() * sense);
1046 let lean = turn.apply_vector(outwards[m]).dot(pull.vector());
1047 if lean > best {
1048 best = lean;
1049 leaning = sense;
1050 }
1051 }
1052 let theta = angle * leaning;
1053
1054 let mut rulings: Vec<Vector> = Vec::with_capacity(hinges.len());
1056 for (hinge, tangent) in hinges.iter().zip(&tangents) {
1057 let turn = Transform::rotation(
1058 ogeom_math::Axis::new(*hinge, Direction::new(*tangent, tol)?),
1059 theta,
1060 );
1061 rulings.push(turn.apply_vector(pull.vector()));
1062 }
1063 let (mut s_lo, mut s_hi) = (f64::INFINITY, f64::NEG_INFINITY);
1070 let (mut h_lo, mut h_hi) = (f64::INFINITY, f64::NEG_INFINITY);
1071 for h in &hinges {
1072 let s = h.to_vector().dot(pull.vector());
1073 h_lo = h_lo.min(s);
1074 h_hi = h_hi.max(s);
1075 }
1076 for p in &mesh.positions {
1077 let s = p.to_vector().dot(pull.vector());
1078 s_lo = s_lo.min(s - h_hi);
1079 s_hi = s_hi.max(s - h_lo);
1080 }
1081 let grow = extent.mul_add(0.5, 1.0) * angle.abs().tan() + tol.confusion();
1082 let (s_lo, s_hi) = (s_lo - grow, s_hi + grow);
1083 let mut lead = 0;
1085 if !closed {
1086 let extend = |hinges: &mut Vec<Point>, rulings: &mut Vec<Vector>, front: bool| -> usize {
1089 let (i0, i1) = if front {
1090 (0, 1)
1091 } else {
1092 (hinges.len() - 1, hinges.len() - 2)
1093 };
1094 let d1 = hinges[i0] - hinges[i1];
1095 let steps = (grow / d1.magnitude().max(tol.confusion())).ceil().max(2.0);
1096 #[allow(clippy::cast_possible_truncation, clippy::cast_sign_loss)]
1097 let steps = (steps as usize).min(16);
1098 for k in 1..=steps {
1099 #[allow(clippy::cast_precision_loss)]
1100 let station = (hinges[i0] + d1 * k as f64, rulings[i0]);
1101 if front {
1102 hinges.insert(0, station.0);
1103 rulings.insert(0, station.1);
1104 } else {
1105 hinges.push(station.0);
1106 rulings.push(station.1);
1107 }
1108 }
1109 steps
1110 };
1111 lead = extend(&mut hinges, &mut rulings, true);
1112 extend(&mut hinges, &mut rulings, false);
1113 } else {
1114 hinges.push(hinges[0]);
1115 rulings.push(rulings[0]);
1116 }
1117 for edge in [s_lo, s_hi] {
1118 for i in 0..hinges.len() - 1 {
1119 let step = (hinges[i + 1] + rulings[i + 1] * edge) - (hinges[i] + rulings[i] * edge);
1120 if step.dot(hinges[i + 1] - hinges[i]) <= 0.0 {
1121 ogeom_bail!(
1122 Construction,
1123 "the draft folds the wall onto itself inside the drafted \
1124 window; refused; see docs/CAPABILITIES.md, offset.draft"
1125 );
1126 }
1127 }
1128 }
1129 let params: Vec<f64> = {
1144 let mut out = Vec::with_capacity(hinges.len());
1145 let mut total = 0.0;
1146 out.push(0.0);
1147 for pair in hinges.windows(2) {
1148 total += pair[0].distance(pair[1]);
1149 out.push(total);
1150 }
1151 if total > 0.0 {
1152 for t in &mut out {
1153 *t /= total;
1154 }
1155 }
1156 if let Some(last) = out.last_mut() {
1157 *last = 1.0;
1158 }
1159 out
1160 };
1161 let stations = feet.len();
1166 let wrap = |x: f64, lo: f64, period: Option<f64>| -> f64 {
1167 period.map_or(x, |p| lo + (x - lo).rem_euclid(p))
1168 };
1169 let exact_at = |at: f64| -> OgeomResult<(Point, Vector)> {
1170 let i = params
1171 .partition_point(|p| *p <= at)
1172 .saturating_sub(1)
1173 .min(params.len() - 2);
1174 let f = (at - params[i]) / (params[i + 1] - params[i]);
1175 if f <= 0.0 {
1176 return Ok((hinges[i], rulings[i]));
1177 }
1178 let (k0, k1) = (i.wrapping_sub(lead), (i + 1).wrapping_sub(lead));
1179 if k0 >= stations || k1 > stations || (!closed && k1 == stations) {
1180 return Ok((hinges[i] + (hinges[i + 1] - hinges[i]) * f, rulings[i]));
1181 }
1182 let (a, b) = (feet[k0], feet[k1 % stations]);
1183 let (mut u, mut v) = (
1184 wrap(a.0 + short(a.0, b.0, period.0) * f, ua, period.0),
1185 wrap(a.1 + short(a.1, b.1, period.1) * f, va, period.1),
1186 );
1187 let (from, chord) = (hinges[i], hinges[i + 1] - hinges[i]);
1192 for _ in 0..8 {
1193 let q = surface.point_at(u, v, tol)?;
1194 let g = (
1195 neutral.signed_distance_to(q),
1196 (q - from).dot(chord) - f * chord.dot(chord),
1197 );
1198 if g.0.abs() <= tol.confusion() * 1e-2
1199 && g.1.abs() <= tol.confusion() * 1e-2 * chord.magnitude()
1200 {
1201 break;
1202 }
1203 let (du, dv) = surface.d1_at(u, v, tol)?;
1204 let (a11, a12, a21, a22) = (n.dot(du), n.dot(dv), chord.dot(du), chord.dot(dv));
1205 let det = a11 * a22 - a12 * a21;
1206 if det.abs() <= f64::MIN_POSITIVE {
1207 break;
1208 }
1209 u = wrap(u - (g.0 * a22 - a12 * g.1) / det, ua, period.0);
1210 v = wrap(v - (a11 * g.1 - a21 * g.0) / det, va, period.1);
1211 }
1212 let (_, p, t, _) = on_hinge((u, v), false)?;
1213 let t = if t.dot(tangents[k0]) < 0.0 { -t } else { t };
1214 let turn = Transform::rotation(ogeom_math::Axis::new(p, Direction::new(t, tol)?), theta);
1215 Ok((p, turn.apply_vector(pull.vector())))
1216 };
1217 const SEAM_BLEND: f64 = 1.0 / 32.0;
1229 let seam = if closed {
1230 let leave = exact_at(params[1] * 1e-9)?.1;
1231 let back = exact_at(1.0 - (1.0 - params[params.len() - 2]) * 1e-9)?.1;
1232 let mean = leave + back;
1233 Some((leave, back, mean / mean.magnitude()))
1234 } else {
1235 None
1236 };
1237 let ruled_at = |at: f64| -> OgeomResult<(Point, Vector)> {
1238 let Some((leave, back, mean)) = seam else {
1239 return exact_at(at);
1240 };
1241 if at <= 0.0 || at >= 1.0 {
1242 return Ok((hinges[0], mean));
1243 }
1244 let (p, r) = exact_at(at)?;
1245 let left = |x: f64| {
1248 let x = x.clamp(0.0, 1.0);
1249 1.0 - x * x * 2.0f64.mul_add(-x, 3.0)
1250 };
1251 let r = if at < SEAM_BLEND {
1252 r + (mean - leave) * left(at / SEAM_BLEND)
1253 } else if at > 1.0 - SEAM_BLEND {
1254 r + (mean - back) * left((1.0 - at) / SEAM_BLEND)
1255 } else {
1256 return Ok((p, r));
1257 };
1258 Ok((p, r / r.magnitude()))
1259 };
1260 let fit_target = (tol.confusion() * 1e3).max(1e-4);
1264 let rows = ogeom_geom::fit::fit_surface_sampled(
1265 |at, side| {
1266 let (h, r) = ruled_at(at)?;
1267 Ok(h + r * if side < 0.5 { s_lo } else { s_hi })
1268 },
1269 &ogeom_geom::fit::Sampling {
1270 us: params.clone(),
1271 vs: vec![0.0, 1.0],
1272 between: (true, false),
1273 closed_v: false,
1274 most: 1024,
1275 },
1276 3,
1277 fit_target,
1278 tol,
1279 )?;
1280 if !rows.met {
1281 ogeom_bail!(
1282 NotDone,
1283 "the drafted wall's border rows reached {} against a target of {fit_target}",
1284 rows.error
1285 );
1286 }
1287 let fitted = rows.curve;
1288 let (k, l) = (fitted.grid().u_count(), fitted.grid().v_count());
1289 if l != 2 || fitted.v_knots().degree() != 1 {
1290 ogeom_bail!(
1291 Construction,
1292 "the drafted wall's border rows did not fit as a ruled surface"
1293 );
1294 }
1295 let corners: Vec<Point> = fitted.grid().points().iter().map(|w| w.point()).collect();
1296 let (lc, hc): (Vec<Point>, Vec<Point>) =
1297 (0..k).map(|i| (corners[i * 2], corners[i * 2 + 1])).unzip();
1298 let mid = {
1303 let (ud, _) = fitted.domain();
1304 f64::midpoint(ud.0, ud.1)
1305 };
1306 let (low_mid, high_mid) = (
1307 fitted.point_at(mid, 0.0, tol)?,
1308 fitted.point_at(mid, 1.0, tol)?,
1309 );
1310 let hinge_mid = low_mid + (high_mid - low_mid) * (-s_lo / (s_hi - s_lo));
1311 let across = fitted.d1_at(mid, 0.0, tol)?.0;
1312 let old_normal = {
1313 let foot = ogeom_algo::project_on_surface(surface, hinge_mid, 16, tol)?;
1314 let (u, v) = foot.parameters;
1315 surface.normal_at(u, v, tol)?.vector()
1316 };
1317 let agrees = across.cross(high_mid - low_mid).dot(old_normal) >= 0.0;
1320 let (first, second, v_range) = if agrees {
1321 (&lc, &hc, (s_lo, s_hi))
1322 } else {
1323 (&hc, &lc, (-s_hi, -s_lo))
1324 };
1325 let mut net: Vec<Point> = Vec::with_capacity(k * 2);
1326 for (a, b) in first.iter().zip(second) {
1327 net.push(*a);
1328 net.push(*b);
1329 }
1330 let grid = ogeom_math::ControlGrid::new(net, k, 2)?;
1331 let u_knots = if closed {
1335 fitted.u_knots().reparameterized(ua, ub)?
1336 } else {
1337 fitted.u_knots().clone()
1338 };
1339 let v_knots =
1340 ogeom_math::KnotVector::clamped_uniform(1, 2)?.reparameterized(v_range.0, v_range.1)?;
1341 Ok(ogeom_geom::BSplineSurface::new(u_knots, v_knots, &grid, tol)?.into())
1342}
1343
1344fn outward_sign(
1353 model: &Model,
1354 solid: &Shape,
1355 face: &Shape,
1356 surface: &SurfaceGeometry,
1357 tol: Tolerances,
1358) -> OgeomResult<f64> {
1359 use ogeom_algo::Containment;
1360 let mesh = ogeom_mesh::triangulate_face(model, face, ogeom_mesh::Deflection::default(), tol)?;
1364 let mut at = None;
1365 let mut largest = 0.0_f64;
1366 for t in &mesh.triangles {
1367 let [a, b, c] = [
1368 mesh.positions[t[0] as usize],
1369 mesh.positions[t[1] as usize],
1370 mesh.positions[t[2] as usize],
1371 ];
1372 let area = (b - a).cross(c - a).magnitude();
1373 if area > largest {
1374 largest = area;
1375 let params = [
1376 mesh.parameters[t[0] as usize],
1377 mesh.parameters[t[1] as usize],
1378 mesh.parameters[t[2] as usize],
1379 ];
1380 at = Some((
1381 (params[0].0 + params[1].0 + params[2].0) / 3.0,
1382 (params[0].1 + params[1].1 + params[2].1) / 3.0,
1383 ));
1384 }
1385 }
1386 let Some((um, vm)) = at else {
1387 ogeom_bail!(Construction, "the drafted face has no interior to probe");
1388 };
1389 let p = surface.point_at(um, vm, tol)?;
1390 let (du, dv) = surface.d1_at(um, vm, tol)?;
1391 let n = du.cross(dv);
1392 let m = n.magnitude();
1393 if m <= tol.confusion() {
1394 ogeom_bail!(Construction, "the face has no normal at its midpoint");
1395 }
1396 let n = n / m;
1397 let scale = largest.sqrt().max(tol.confusion() * 1e3);
1398 for eps_scale in [1e-3, 1e-2, 5e-2] {
1399 let eps = scale * eps_scale;
1400 let deflection = ogeom_mesh::Deflection {
1401 chord: (eps * 0.1).max(1e-4),
1402 ..ogeom_mesh::Deflection::default()
1403 };
1404 let ahead = ogeom_algo::classify_in_solid(model, solid, p + n * eps, deflection, tol)?;
1405 let behind = ogeom_algo::classify_in_solid(model, solid, p - n * eps, deflection, tol)?;
1406 match (ahead, behind) {
1407 (Containment::Out, Containment::In) => return Ok(1.0),
1408 (Containment::In, Containment::Out) => return Ok(-1.0),
1409 _ => {}
1410 }
1411 }
1412 ogeom_bail!(
1413 Construction,
1414 "cannot read which side of the drafted face holds material; the wall is thinner than the probe can resolve"
1415 )
1416}
1417
1418fn meet(a: Plane, b: Plane, along: Vector, tol: Tolerances) -> OgeomResult<Point> {
1420 let rows = [a.normal().vector(), b.normal().vector(), along];
1421 let rhs = [
1422 rows[0].dot(a.origin().to_vector()),
1423 rows[1].dot(b.origin().to_vector()),
1424 along.dot(Point::midpoint(a.origin(), b.origin()).to_vector()),
1425 ];
1426 let det = rows[0].dot(rows[1].cross(rows[2]));
1427 if det.abs() <= tol.confusion() {
1428 ogeom_bail!(Construction, "the two planes do not meet in a line");
1429 }
1430 Ok(Point::ORIGIN
1431 + (rows[1].cross(rows[2]) * rhs[0]
1432 + rows[2].cross(rows[0]) * rhs[1]
1433 + rows[0].cross(rows[1]) * rhs[2])
1434 / det)
1435}