1use super::*;
135use crate::offset_retrim::plane_of_surface;
136use crate::topology::ShellRecord;
137use crate::{
138 boolean_operation, make_plane, BooleanOperation, BooleanOptions,
139
140};
141use crate::curve::parameter_line;
142use crate::{KernelRefusal, KernelStage, OrRefuse};
143
144const OP: &str = "move_faces(re-cut)";
145
146#[derive(Clone, Copy, Debug)]
148struct HalfSpace {
149 normal: Vec3,
150 offset: f64,
151}
152
153
154pub fn recut_moved_planes(
161 solid: &BrepSolid,
162 moved: &[u64],
163 translation: Vec3,
164) -> Result<BrepSolid, KernelRefusal> {
165 recut_moved_planes_rigid(solid, moved, None, translation)
166}
167
168struct Motion {
172 rotation: Option<AffineTransform>,
173 translation: Vec3,
174}
175
176impl Motion {
177 fn point(&self, point: Vec3) -> Vec3 {
178 match &self.rotation {
179 None => point.add(self.translation),
180 Some(rotation) => rotation.point(point).add(self.translation),
181 }
182 }
183
184 fn plane(&self, normal: Vec3, offset: f64) -> (Vec3, f64) {
188 match &self.rotation {
189 None => (normal, offset + self.translation.dot(normal)),
190 Some(rotation) => {
191 let turned = super::face_rotate::rotate_direction(rotation, normal);
192 (turned, turned.dot(self.point(normal.scale(offset))))
193 }
194 }
195 }
196}
197
198pub fn recut_moved_planes_rigid(
222 solid: &BrepSolid,
223 moved: &[u64],
224 rotation: Option<(Vec3, Vec3, f64)>,
225 translation: Vec3,
226) -> Result<BrepSolid, KernelRefusal> {
227 let scale = solid_model_scale(solid);
228 let tolerance = (scale * 1e-7).max(1e-9);
229 let moved_set: HashSet<u64> = moved.iter().copied().collect();
230 if moved_set.is_empty() {
231 return Err(KernelRefusal::input(
232 KernelStage::Collect,
233 "selection",
234 format!("{OP}: no face was selected"),
235 ));
236 }
237 let motion = Motion {
238 rotation: match rotation {
239 None => None,
240 Some((point, axis, angle)) => {
241 let axis = axis.normalized().map_err(|_| {
242 KernelRefusal::input(
243 KernelStage::Collect,
244 "rotation_axis",
245 format!("{OP}: the rotation axis has zero length"),
246 )
247 })?;
248 Some(
249 super::face_rotate::rotation_about(point, axis, angle)
250 .or_input(KernelStage::Collect, "rotation")?,
251 )
252 }
253 },
254 translation,
255 };
256
257 let mut carriers: Vec<(u64, Vec3, f64)> = Vec::new();
261 for &id in &moved_set {
262 let face = face_by_id(solid, id)?;
263 let plane = plane_of_surface(&face.surface, tolerance.max(1e-9) * 10.0, OP).map_err(
264 |error| {
265 KernelRefusal::unsupported(
266 KernelStage::Classify,
267 "curved_carrier",
268 format!(
269 "{error} — this road sweeps a PLANE's slab, so face {id}{} cannot take it",
270 named(face)
271 ),
272 )
273 },
274 )?;
275 let normal = if face.same_sense {
276 plane.normal
277 } else {
278 plane.normal.scale(-1.0)
279 };
280 carriers.push((id, normal, normal.dot(plane.origin)));
281 }
282 carriers.sort_by_key(|entry| entry.0);
284
285 let bounds = solid_bounds(solid)?;
290 let bounds = if motion.rotation.is_none() {
291 bounds
292 } else {
293 solid.vertices.iter().fold(bounds, |(low, high), vertex| {
294 let point = motion.point(vertex.point);
295 (
296 Vec3::new(low.x.min(point.x), low.y.min(point.y), low.z.min(point.z)),
297 Vec3::new(high.x.max(point.x), high.y.max(point.y), high.z.max(point.z)),
298 )
299 })
300 };
301 let margin = translation.length() + scale.max(1.0);
302 let clip = box_half_spaces(bounds, margin);
303
304 let motion_floor = (scale * 1e-9).max(1e-12);
305 let edge_band = (scale * 1e-7).max(1e-9);
306 let mut removing: Vec<BrepSolid> = Vec::new();
307 let mut adding: Vec<BrepSolid> = Vec::new();
308
309 for &(id, normal, offset) in &carriers {
310 let face = face_by_id(solid, id)?;
311 let (moved_normal, moved_offset) = motion.plane(normal, offset);
312 let facing = moved_normal.dot(normal);
313 if motion.rotation.is_some() && facing <= 1e-6 {
314 return Err(KernelRefusal::ill_posed(
315 KernelStage::Classify,
316 "turned_past_right_angle",
317 format!(
318 "{OP}: the motion turns face {id}{} through a right angle or more, so its \
319 new carrier does not face the way the old one did and there is no wedge \
320 between them for a tool to fill; refusing",
321 named(face)
322 ),
323 ));
324 }
325 let unturned = moved_normal.sub(normal).length() <= 1e-12;
326 let translation_delta = translation.dot(normal);
329 if motion.rotation.is_none() && translation_delta.abs() <= motion_floor {
330 continue;
332 }
333
334 let mut neighbour_planes: HashMap<u64, ((Vec3, f64), (Vec3, f64))> = HashMap::default();
337 let mut tangent_bands: HashSet<u64> = HashSet::default();
339 let polygon = trim_polygon(solid, face, normal, edge_band)?;
340 for source in &polygon.sources {
341 let &EdgeSource::Neighbour(neighbour) = source else { continue };
342 if neighbour_planes.contains_key(&neighbour) || tangent_bands.contains(&neighbour) {
343 continue;
344 }
345 let other = face_by_id(solid, neighbour)?;
346 if motion.rotation.is_some()
347 && !moved_set.contains(&neighbour)
348 && plane_of_surface(&other.surface, tolerance.max(1e-9) * 10.0, OP).is_err()
349 && band_is_tangent_along_shared_edges(face, other, normal)
350 {
351 tangent_bands.insert(neighbour);
352 continue;
353 }
354 let plane = plane_of_surface(&other.surface, tolerance.max(1e-9) * 10.0, OP)
355 .map_err(|error| {
356 KernelRefusal::unsupported(
357 KernelStage::Classify,
358 "curved_neighbour",
359 format!(
360 "{error} — face {id}{} is bounded by face {neighbour}{}, and this road \
361 needs every bounding carrier as a half-space",
362 named(face),
363 named(other)
364 ),
365 )
366 })?;
367 let outward = if other.same_sense {
368 plane.normal
369 } else {
370 plane.normal.scale(-1.0)
371 };
372 let home = (outward, outward.dot(plane.origin));
373 let left = if moved_set.contains(&neighbour) {
374 motion.plane(home.0, home.1)
375 } else {
376 home
377 };
378 neighbour_planes.insert(neighbour, (home, left));
379 }
380
381 let hinge = if unturned {
384 None
385 } else {
386 polygon.frame.hinge(normal, offset, moved_normal, moved_offset)
387 };
388 let cells = convex_cells(&polygon, hinge, edge_band, scale).map_err(|error| {
389 KernelRefusal::internal(
390 KernelStage::Classify,
391 "cell_split",
392 format!("{OP}: face {id}{}'s trim could not be cut into convex cells — {error}", named(face)),
393 )
394 })?;
395
396 let surround = if tangent_bands.is_empty() {
397 None
398 } else {
399 Some(tangent_surround(solid, face, normal, offset, &tangent_bands, &moved_set, tolerance)?)
400 };
401
402 for cell in &cells {
403 let seat = polygon.frame.lift(cell.centroid());
404 let delta = if motion.rotation.is_none() {
409 translation_delta
410 } else {
411 (moved_offset - moved_normal.dot(seat)) / facing
412 };
413 if delta.abs() <= motion_floor {
414 if unturned {
415 continue;
416 }
417 return Err(KernelRefusal::ill_posed(
418 KernelStage::Classify,
419 "turned_through_trim",
420 format!(
421 "{OP}: the motion turns face {id}{} about a line through its own trim, and a \
422 cell of the trim still sits on that line after the hinge split, so which way \
423 it moves is not decided by the geometry; refusing rather than choosing",
424 named(face)
425 ),
426 ));
427 }
428 let adds = delta > 0.0;
429
430 let (inner, outer) = if adds {
435 ((moved_normal, moved_offset), (normal, offset))
436 } else {
437 ((normal, offset), (moved_normal, moved_offset))
438 };
439 let wedge = |inner: (Vec3, f64), outer: (Vec3, f64)| {
440 [
441 HalfSpace { normal: inner.0, offset: inner.1 },
442 HalfSpace { normal: outer.0.scale(-1.0), offset: -outer.1 },
443 ]
444 };
445 let touches_band = cell.sources.iter().any(|source| {
448 matches!(source, EdgeSource::Neighbour(neighbour) if tangent_bands.contains(neighbour))
449 });
450 let (inner, outer) = match (&surround, adds && touches_band) {
451 (Some(surround), true) => (inner, (outer.0, outer.1 - surround.inner_wall_depth)),
452 _ => (inner, outer),
453 };
454 let mut spaces = clip.clone();
455 if let (Some(surround), true) = (&surround, touches_band) {
456 for plane in &surround.planes {
457 push_unique(
458 &mut spaces,
459 HalfSpace { normal: plane.0, offset: plane.1 },
460 tolerance,
461 );
462 }
463 }
464 if motion.rotation.is_none() {
465 spaces.push(HalfSpace {
468 normal,
469 offset: offset.max(offset + delta),
470 });
471 spaces.push(HalfSpace {
472 normal: normal.scale(-1.0),
473 offset: -offset.min(offset + delta),
474 });
475 } else {
476 spaces.extend(wedge(inner, outer));
477 }
478
479 let probe = seat.add(normal.scale(delta * 0.5));
483
484 let mut at_home: Vec<HalfSpace> = Vec::new();
487 let mut sides: Vec<HalfSpace> = Vec::new();
489
490 for (edge, source) in cell.sources.iter().enumerate() {
491 if matches!(source, EdgeSource::Neighbour(neighbour) if tangent_bands.contains(neighbour)) {
495 let plane = polygon.frame.edge_plane(cell, edge);
496 let home_sign = if plane.0.dot(seat) - plane.1 > 0.0 { -1.0 } else { 1.0 };
497 push_unique(
498 &mut at_home,
499 HalfSpace { normal: plane.0.scale(home_sign), offset: plane.1 * home_sign },
500 tolerance,
501 );
502 continue;
503 }
504 let (home, bound) = match source {
505 EdgeSource::Neighbour(neighbour) => {
506 let &(home, left) = neighbour_planes
507 .get(neighbour)
508 .expect("every neighbour's plane was read above");
509 let material_side = home.0.dot(seat) - home.1 < 0.0;
523 let bound = if material_side == adds { left } else { home };
524 (home, bound)
525 }
526 EdgeSource::Hinge => {
527 let plane = polygon.frame.edge_plane(cell, edge);
528 (plane, plane)
529 }
530 };
531 let reach = bound.0.dot(probe) - bound.1;
532 if reach.abs() <= tolerance {
533 return Err(KernelRefusal::ill_posed(
534 KernelStage::Classify,
535 "probe_on_carrier",
536 format!(
537 "{OP}: the probe point for face {id} lies ON the carrier of one of its \
538 bounds ({reach:.3e} from it), so which side of that carrier the tool \
539 is on is not decided by the geometry; refusing rather than choosing"
540 ),
541 ));
542 }
543 let sign = if reach > 0.0 { -1.0 } else { 1.0 };
544 let bound = HalfSpace {
545 normal: bound.0.scale(sign),
546 offset: bound.1 * sign,
547 };
548 let home_sign = if home.0.dot(seat) - home.1 > 0.0 { -1.0 } else { 1.0 };
551 push_unique(&mut spaces, bound, tolerance);
552 push_unique(&mut sides, bound, tolerance);
553 push_unique(
554 &mut at_home,
555 HalfSpace {
556 normal: home.0.scale(home_sign),
557 offset: home.1 * home_sign,
558 },
559 tolerance,
560 );
561 }
562 cell_is_its_bounds(face, cell, &polygon.frame, &at_home, scale)?;
563 if motion.rotation.is_some() {
564 let (home_inner, home_outer) = if adds {
565 ((moved_normal, moved_offset), (normal, offset))
566 } else {
567 ((normal, offset), (moved_normal, moved_offset))
568 };
569 one_way_only(id, &clip, &sides, wedge(home_outer, home_inner), tolerance, || named(face))?;
570 }
571
572 let tool = convex_solid(&spaces, probe, tolerance).map_err(|error| {
573 KernelRefusal::internal(
574 KernelStage::Fragment,
575 "tool_polytope",
576 format!("{OP}: the tool for face {id} could not be built — {error}"),
577 )
578 })?;
579 if adds {
580 adding.push(tool);
581 } else {
582 removing.push(tool);
583 }
584 }
585 }
586 if removing.is_empty() && adding.is_empty() {
587 return Err(KernelRefusal::input(
588 KernelStage::Classify,
589 "in_plane_motion",
590 format!(
591 "{OP}: every selected face slides inside its own carrier (the motion has no \
592 component along any of their normals), so there is no material for a re-cut to \
593 move"
594 ),
595 ));
596 }
597
598 let unite = |mut tools: Vec<BrepSolid>| -> Result<Option<BrepSolid>, KernelRefusal> {
606 if tools.is_empty() {
607 return Ok(None);
608 }
609 let mut tool = tools.remove(0);
610 for next in tools {
611 tool = boolean_operation(&tool, &next, BooleanOperation::Union, &BooleanOptions::default())
614 .map_err(|error| {
615 error.with_message(|message| {
616 format!("{OP}: the tools could not be unioned into one — {message}")
617 })
618 })?;
619 }
620 Ok(Some(tool))
621 };
622 let mut cut = solid.clone();
623 if let Some(tool) = unite(removing)? {
624 cut = boolean_operation(&cut, &tool, BooleanOperation::Subtract, &BooleanOptions::default())
625 .map_err(|error| {
626 error.with_message(|message| format!("{OP}: the re-cut's subtract refused — {message}"))
627 })?;
628 }
629 if let Some(tool) = unite(adding)? {
630 cut = boolean_operation(&cut, &tool, BooleanOperation::Union, &BooleanOptions::default())
631 .map_err(|error| {
632 error.with_message(|message| format!("{OP}: the re-cut's union refused — {message}"))
633 })?;
634 }
635 Ok(restamp(cut, solid, &moved_set, &motion, tolerance))
636}
637
638pub fn reblend_lost_fillets(original: &BrepSolid, moved: &[u64], mut body: BrepSolid) -> (BrepSolid, Vec<String>) {
646 let tolerance = (solid_model_scale(original) * 1e-6).max(1e-9);
647 let moved_set: HashSet<u64> = moved.iter().copied().collect();
648 let faces: HashMap<u64, &FaceRecord> =
649 original.shells.iter().flat_map(|shell| &shell.faces).map(|face| (face.id, face)).collect();
650 let mut owners: HashMap<u64, Vec<u64>> = HashMap::default();
651 for face in faces.values() {
652 for coedge in face.loops.iter().flat_map(|wire| &wire.coedges) {
653 owners.entry(coedge.edge_id).or_default().push(face.id);
654 }
655 }
656 let across = |face: &FaceRecord| -> Vec<u64> {
657 let mut found: Vec<u64> = face
658 .loops
659 .iter()
660 .flat_map(|wire| &wire.coedges)
661 .flat_map(|coedge| owners.get(&coedge.edge_id).into_iter().flatten().copied())
662 .filter(|id| *id != face.id)
663 .collect();
664 found.sort_unstable();
665 found.dedup();
666 found
667 };
668 let mut notes = Vec::new();
669 let mut ordered = moved.to_vec();
670 ordered.sort_unstable();
671 for id in ordered {
672 let Some(&face) = faces.get(&id) else { continue };
673 let Ok(plane) = plane_of_surface(&face.surface, tolerance * 10.0, OP) else { continue };
674 for band_id in across(face) {
675 let band = faces[&band_id];
676 let Some(band_name) = band.name.clone() else { continue };
677 if moved_set.contains(&band_id)
678 || plane_of_surface(&band.surface, tolerance * 10.0, OP).is_ok()
679 || !band_is_tangent_along_shared_edges(face, band, plane.normal)
680 || body.shells.iter().flat_map(|shell| &shell.faces).any(|f| f.name.as_deref() == Some(&band_name))
681 {
682 continue;
683 }
684 let radius = super::face_move::ruled_revolution_carrier(&band.surface)
687 .filter(|(_, rho0, rho1, _)| (rho0 - rho1).abs() <= tolerance)
688 .map(|(_, rho0, _, _)| rho0);
689 let Some(rho0) = radius else {
690 notes.push(format!(
691 "transformFace: fillet `{band_name}` was lost: the moved face no longer reaches \
692 it, so it was filled in, and only a constant-radius fillet can be re-made"
693 ));
694 continue;
695 };
696 let remade = across(band)
697 .into_iter()
698 .filter(|wall| *wall != id)
699 .find_map(|wall| {
700 let wall = faces[&wall];
701 let wall_plane = plane_of_surface(&wall.surface, tolerance * 10.0, OP).ok()?;
702 band_is_tangent_along_shared_edges(wall, band, wall_plane.normal).then_some(wall)
703 })
704 .ok_or_else(|| "the wall beyond it could not be read".to_string())
705 .and_then(|wall| sharp_edge_between(&body, face, wall))
706 .and_then(|edge| {
707 let rebuilt = crate::fillet_edge(&body, edge, rho0, Some(&band_name)).map_err(|e| e.to_string())?;
708 let issues = rebuilt.validate();
709 if !issues.is_empty() {
710 return Err(format!("the new fillet is invalid: {issues:?}"));
711 }
712 crate::accept_sound(rebuilt, "transformFace").map_err(|e| e.to_string())
715 });
716 match remade {
717 Ok(rebuilt) => {
718 body = rebuilt;
719 for beside in across(band) {
722 let beside = faces[&beside];
723 let Some(beside_name) = beside.name.as_deref() else { continue };
724 if !moved_set.contains(&beside.id)
725 && plane_of_surface(&beside.surface, tolerance * 10.0, OP).is_err()
726 && !body.shells.iter().flat_map(|shell| &shell.faces).any(|f| f.name.as_deref() == Some(beside_name))
727 {
728 notes.push(format!(
729 "transformFace: blend `{beside_name}`, which ran round fillet \
730 `{band_name}`, was not re-made: the fillet was rebuilt at its radius \
731 and the faces either side of that blend now meet over it"
732 ));
733 }
734 }
735 }
736 Err(reason) => notes.push(format!(
737 "transformFace: fillet `{band_name}` (radius {rho0}) was lost: the moved face no \
738 longer reaches it, so it was filled in, and it could not be re-made ({reason})"
739 )),
740 }
741 }
742 }
743 (body, notes)
744}
745
746fn sharp_edge_between(body: &BrepSolid, first: &FaceRecord, second: &FaceRecord) -> Result<u64, String> {
748 let edges_of = |face: &FaceRecord| -> Result<HashSet<u64>, String> {
749 let name = face.name.as_deref().ok_or_else(|| format!("face {} has no name to find it by", face.id))?;
750 Ok(body
751 .shells
752 .iter()
753 .flat_map(|shell| &shell.faces)
754 .filter(|f| f.name.as_deref() == Some(name))
755 .flat_map(|f| f.loops.iter().flat_map(|wire| &wire.coedges).map(|coedge| coedge.edge_id))
756 .collect())
757 };
758 let (a, b) = (edges_of(first)?, edges_of(second)?);
759 let mut shared = a.intersection(&b).copied();
760 match (shared.next(), shared.next()) {
761 (Some(edge), None) => Ok(edge),
762 _ => Err("no single sharp edge joins the turned face to the wall beyond it".to_string()),
763 }
764}
765
766struct TangentSurround {
768 planes: Vec<(Vec3, f64)>,
770 inner_wall_depth: f64,
772}
773
774fn band_is_tangent_along_shared_edges(
776 face: &FaceRecord,
777 band: &FaceRecord,
778 normal: Vec3,
779) -> bool {
780 let own: HashSet<u64> =
781 face.loops.iter().flat_map(|wire| &wire.coedges).map(|coedge| coedge.edge_id).collect();
782 let mut shared = 0usize;
783 for coedge in band.loops.iter().flat_map(|wire| &wire.coedges) {
784 if !own.contains(&coedge.edge_id) {
785 continue;
786 }
787 let Ok([q0, q1]) = coedge.pcurve.domain() else {
788 return false;
789 };
790 let Ok(uv) = coedge.pcurve.evaluate(0.5 * (q0 + q1)) else {
791 return false;
792 };
793 let Ok(band_normal) = band.surface.normal(uv.x, uv.y) else {
794 return false;
795 };
796 if band_normal.cross(normal).length() > 1e-4 {
797 return false;
798 }
799 shared += 1;
800 }
801 shared > 0
802}
803
804fn tangent_surround(
813 solid: &BrepSolid,
814 face: &FaceRecord,
815 normal: Vec3,
816 offset: f64,
817 bands: &HashSet<u64>,
818 moved: &HashSet<u64>,
819 tolerance: f64,
820) -> Result<TangentSurround, KernelRefusal> {
821 let mut owners: HashMap<u64, Vec<u64>> = HashMap::default();
822 for other in solid.shells.iter().flat_map(|shell| &shell.faces) {
823 for coedge in other.loops.iter().flat_map(|wire| &wire.coedges) {
824 owners.entry(coedge.edge_id).or_default().push(other.id);
825 }
826 }
827 let edges: HashMap<u64, &EdgeRecord> = solid.edges.iter().map(|edge| (edge.id, edge)).collect();
828 let points: HashMap<u64, Vec3> =
829 solid.vertices.iter().map(|vertex| (vertex.id, vertex.point)).collect();
830 let plane_band = tolerance.max(1e-9) * 10.0;
831
832 let mut curved: HashSet<u64> = bands.clone();
833 let mut planar: HashSet<u64> = HashSet::default();
834 let mut stack: Vec<u64> = bands.iter().copied().collect();
835 while let Some(id) = stack.pop() {
836 let current = face_by_id(solid, id)?;
837 for coedge in current.loops.iter().flat_map(|wire| &wire.coedges) {
838 for &next in owners.get(&coedge.edge_id).into_iter().flatten() {
839 if next == id || next == face.id || curved.contains(&next) || planar.contains(&next) {
840 continue;
841 }
842 if moved.contains(&next) {
843 return Err(KernelRefusal::unsupported(
844 KernelStage::Classify,
845 "curved_neighbour",
846 format!(
847 "{OP}: face {next} is also selected and borders the fillet beside face \
848 {}, so what bounds the tool past that fillet is itself moving; refusing",
849 face.id
850 ),
851 ));
852 }
853 let other = face_by_id(solid, next)?;
854 if plane_of_surface(&other.surface, plane_band, OP).is_ok() {
855 planar.insert(next);
856 } else {
857 curved.insert(next);
858 stack.push(next);
859 }
860 }
861 }
862 }
863
864 let vertices_of = |id: u64| -> Result<Vec<Vec3>, KernelRefusal> {
865 let mut found = Vec::new();
866 for coedge in face_by_id(solid, id)?.loops.iter().flat_map(|wire| &wire.coedges) {
867 if let Some(edge) = edges.get(&coedge.edge_id) {
868 found.extend(points.get(&edge.start_vertex_id).copied());
869 found.extend(points.get(&edge.end_vertex_id).copied());
870 }
871 }
872 Ok(found)
873 };
874 let mut region: Vec<Vec3> = Vec::new();
875 let mut depth = 0.0f64;
876 for &id in &curved {
877 for point in vertices_of(id)? {
878 depth = depth.max(offset - normal.dot(point));
879 region.push(point);
880 }
881 }
882 let mut ordered: Vec<u64> = planar.iter().copied().collect();
883 ordered.sort_unstable();
884 for &id in &ordered {
885 region.extend(vertices_of(id)?);
886 }
887
888 let mut planes = Vec::new();
889 for &id in &ordered {
890 let other = face_by_id(solid, id)?;
891 let plane = plane_of_surface(&other.surface, plane_band, OP)
892 .map_err(|error| KernelRefusal::unsupported(KernelStage::Classify, "curved_neighbour", error.to_string()))?;
893 let outward = if other.same_sense { plane.normal } else { plane.normal.scale(-1.0) };
894 let outward_offset = outward.dot(plane.origin);
895 if let Some(point) = region.iter().find(|point| outward.dot(**point) - outward_offset > plane_band) {
896 return Err(KernelRefusal::unsupported(
897 KernelStage::Classify,
898 "curved_neighbour",
899 format!(
900 "{OP}: the part is not convex around the fillet beside face {} (face {id}{} has \
901 material on both sides, at {point:?}), so the wedge a turn adds there has no \
902 bound that keeps it out of the part's other features; refusing",
903 face.id,
904 named(other)
905 ),
906 ));
907 }
908 planes.push((outward, outward_offset));
909 }
910 Ok(TangentSurround { planes, inner_wall_depth: depth * 1.5 + tolerance })
911}
912
913fn push_unique(spaces: &mut Vec<HalfSpace>, space: HalfSpace, tolerance: f64) {
918 let band = tolerance.max(1e-9) * 10.0;
919 if spaces.iter().any(|other| {
920 other.normal.sub(space.normal).length() <= 1e-9 && (other.offset - space.offset).abs() <= band
921 }) {
922 return;
923 }
924 spaces.push(space);
925}
926
927#[derive(Clone, Copy, Debug, PartialEq)]
931enum EdgeSource {
932 Neighbour(u64),
933 Hinge,
934}
935
936#[derive(Clone, Copy, Debug)]
940struct Frame {
941 origin: Vec3,
942 u: Vec3,
943 v: Vec3,
944}
945
946impl Frame {
947 fn new(origin: Vec3, normal: Vec3) -> Result<Self, String> {
948 let u = normal.perpendicular()?.normalized()?;
949 let v = normal.cross(u).normalized()?;
950 Ok(Self { origin, u, v })
951 }
952
953 fn drop(&self, point: Vec3) -> [f64; 2] {
954 let delta = point.sub(self.origin);
955 [delta.dot(self.u), delta.dot(self.v)]
956 }
957
958 fn lift(&self, point: [f64; 2]) -> Vec3 {
959 self.origin.add(self.u.scale(point[0])).add(self.v.scale(point[1]))
960 }
961
962 fn hinge(&self, normal: Vec3, offset: f64, moved: Vec3, moved_offset: f64) -> Option<Line2> {
965 let direction = normal.cross(moved);
966 if direction.length() <= 1e-12 {
967 return None;
968 }
969 let direction = direction.normalized().ok()?;
970 let across = direction.cross(normal).normalized().ok()?;
973 let rate = moved.dot(across);
974 if rate.abs() <= 1e-12 {
975 return None;
976 }
977 let _ = offset;
978 let point = self.origin.add(across.scale((moved_offset - moved.dot(self.origin)) / rate));
979 let [px, py] = self.drop(point);
980 let [qx, qy] = self.drop(point.add(direction));
981 Some(Line2 {
982 point: [px, py],
983 direction: [qx - px, qy - py],
984 })
985 }
986
987 fn edge_plane(&self, cell: &Cell, edge: usize) -> (Vec3, f64) {
991 let a = cell.points[edge];
992 let b = cell.points[(edge + 1) % cell.points.len()];
993 let d = [b[0] - a[0], b[1] - a[1]];
996 let length = (d[0] * d[0] + d[1] * d[1]).sqrt().max(1e-300);
997 let out = [d[1] / length, -d[0] / length];
998 let normal = self.u.scale(out[0]).add(self.v.scale(out[1]));
999 let normal = normal.normalized().unwrap_or(normal);
1000 (normal, normal.dot(self.lift(a)))
1001 }
1002}
1003
1004#[derive(Clone, Copy, Debug)]
1006struct Line2 {
1007 point: [f64; 2],
1008 direction: [f64; 2],
1009}
1010
1011impl Line2 {
1012 fn side(&self, p: [f64; 2]) -> f64 {
1015 let d = self.direction;
1016 let length = (d[0] * d[0] + d[1] * d[1]).sqrt().max(1e-300);
1017 (d[0] * (p[1] - self.point[1]) - d[1] * (p[0] - self.point[0])) / length
1018 }
1019}
1020
1021#[derive(Clone, Debug)]
1025struct Cell {
1026 points: Vec<[f64; 2]>,
1027 sources: Vec<EdgeSource>,
1028}
1029
1030impl Cell {
1031 fn signed_area(&self) -> f64 {
1032 let n = self.points.len();
1033 (0..n)
1034 .map(|i| {
1035 let a = self.points[i];
1036 let b = self.points[(i + 1) % n];
1037 crate::polygon::area_term(a, b)
1038 })
1039 .sum::<f64>()
1040 * 0.5
1041 }
1042
1043 fn centroid(&self) -> [f64; 2] {
1045 let n = self.points.len();
1046 let area = self.signed_area();
1047 if area.abs() <= 1e-300 {
1048 let mean = self.points.iter().fold([0.0, 0.0], |s, p| [s[0] + p[0], s[1] + p[1]]);
1049 return [mean[0] / n as f64, mean[1] / n as f64];
1050 }
1051 let mut cx = 0.0;
1052 let mut cy = 0.0;
1053 for i in 0..n {
1054 let a = self.points[i];
1055 let b = self.points[(i + 1) % n];
1056 let w = crate::polygon::area_term(a, b);
1057 cx += (a[0] + b[0]) * w;
1058 cy += (a[1] + b[1]) * w;
1059 }
1060 [cx / (6.0 * area), cy / (6.0 * area)]
1061 }
1062
1063 fn reflex_vertex(&self, band: f64) -> Option<usize> {
1067 let n = self.points.len();
1068 (0..n).find(|&i| {
1069 let p = self.points[(i + n - 1) % n];
1070 let q = self.points[i];
1071 let r = self.points[(i + 1) % n];
1072 let a = [q[0] - p[0], q[1] - p[1]];
1073 let b = [r[0] - q[0], r[1] - q[1]];
1074 let la = (a[0] * a[0] + a[1] * a[1]).sqrt();
1075 let lb = (b[0] * b[0] + b[1] * b[1]).sqrt();
1076 if la <= band || lb <= band {
1077 return false;
1078 }
1079 (a[0] * b[1] - a[1] * b[0]) / la.max(lb) < -band
1082 })
1083 }
1084}
1085
1086struct TrimPolygon {
1088 frame: Frame,
1089 points: Vec<[f64; 2]>,
1090 sources: Vec<EdgeSource>,
1091}
1092
1093fn trim_polygon(
1100 solid: &BrepSolid,
1101 face: &FaceRecord,
1102 normal: Vec3,
1103 band: f64,
1104) -> Result<TrimPolygon, KernelRefusal> {
1105 let id = face.id;
1106 if face.loops.len() != 1 {
1107 return Err(KernelRefusal::unsupported(
1108 KernelStage::Classify,
1109 "island_trim",
1110 format!(
1111 "{OP}: face {id}{} has {} loops, so its trim has an island whose walls bound no cell \
1112 of the outer trim, and this road cuts one convex cell at a time; refusing rather \
1113 than filling the island",
1114 named(face),
1115 face.loops.len()
1116 ),
1117 ));
1118 }
1119 let edges: HashMap<u64, &EdgeRecord> = solid.edges.iter().map(|edge| (edge.id, edge)).collect();
1120 let vertices: HashMap<u64, Vec3> =
1121 solid.vertices.iter().map(|vertex| (vertex.id, vertex.point)).collect();
1122 let mut across: HashMap<u64, u64> = HashMap::default();
1124 for other in solid.shells.iter().flat_map(|shell| &shell.faces) {
1125 if other.id == id {
1126 continue;
1127 }
1128 for coedge in other.loops.iter().flat_map(|wire| &wire.coedges) {
1129 across.entry(coedge.edge_id).or_insert(other.id);
1130 }
1131 }
1132 let mut points3: Vec<Vec3> = Vec::new();
1133 let mut sources: Vec<EdgeSource> = Vec::new();
1134 for coedge in &face.loops[0].coedges {
1135 let edge = edges.get(&coedge.edge_id).ok_or_else(|| {
1136 KernelRefusal::input(
1137 KernelStage::Collect,
1138 "missing_edge",
1139 format!("{OP}: face {id} uses edge {} which the solid does not hold", coedge.edge_id),
1140 )
1141 })?;
1142 if edge.degenerate {
1143 continue;
1144 }
1145 let (start_id, end_id) = if coedge.forward {
1146 (edge.start_vertex_id, edge.end_vertex_id)
1147 } else {
1148 (edge.end_vertex_id, edge.start_vertex_id)
1149 };
1150 let start = *vertices.get(&start_id).ok_or_else(|| {
1151 KernelRefusal::input(
1152 KernelStage::Collect,
1153 "missing_vertex",
1154 format!("{OP}: edge {} starts at vertex {start_id}, which the solid does not hold", edge.id),
1155 )
1156 })?;
1157 let end = *vertices.get(&end_id).ok_or_else(|| {
1158 KernelRefusal::input(
1159 KernelStage::Collect,
1160 "missing_vertex",
1161 format!("{OP}: edge {} ends at vertex {end_id}, which the solid does not hold", edge.id),
1162 )
1163 })?;
1164 let middle = edge
1166 .curve
1167 .evaluate(0.5 * (edge.t0 + edge.t1))
1168 .or_refuse(KernelStage::Classify, "evaluate")?;
1169 let chord_middle = start.add(end).scale(0.5);
1170 if middle.sub(chord_middle).length() > band.max(1e-9) * 100.0 {
1171 return Err(KernelRefusal::unsupported(
1172 KernelStage::Classify,
1173 "curved_edge",
1174 format!(
1175 "{OP}: edge {} of face {id}{} is not straight (its midpoint is {:.3e} off its \
1176 chord), so the trim is not a polygon this road can cut into cells",
1177 edge.id,
1178 named(face),
1179 middle.sub(chord_middle).length()
1180 ),
1181 ));
1182 }
1183 let neighbour = *across.get(&edge.id).ok_or_else(|| {
1184 KernelRefusal::unsupported(
1185 KernelStage::Classify,
1186 "open_edge",
1187 format!(
1188 "{OP}: edge {} of face {id}{} is held by no other face, so nothing bounds the \
1189 tool there",
1190 edge.id,
1191 named(face)
1192 ),
1193 )
1194 })?;
1195 points3.push(start);
1196 sources.push(EdgeSource::Neighbour(neighbour));
1197 }
1198 if points3.len() < 3 {
1199 return Err(KernelRefusal::unsupported(
1200 KernelStage::Classify,
1201 "thin_trim",
1202 format!("{OP}: face {id}{} has fewer than three corners", named(face)),
1203 ));
1204 }
1205 let frame = Frame::new(points3[0], normal).map_err(|error| {
1206 KernelRefusal::internal(KernelStage::Classify, "frame", format!("{OP}: {error}"))
1207 })?;
1208 let mut points: Vec<[f64; 2]> = points3.iter().map(|point| frame.drop(*point)).collect();
1209 let probe = Cell { points: points.clone(), sources: sources.clone() };
1210 if probe.signed_area() < 0.0 {
1211 points.reverse();
1214 sources.reverse();
1215 sources.rotate_left(1);
1216 let _ = &mut points;
1217 }
1218 Ok(TrimPolygon { frame, points, sources })
1219}
1220
1221fn convex_cells(
1227 polygon: &TrimPolygon,
1228 hinge: Option<Line2>,
1229 band: f64,
1230 scale: f64,
1231) -> Result<Vec<Cell>, String> {
1232 let area_floor = (band * scale).max(1e-18);
1233 let mut pending = vec![Cell {
1234 points: polygon.points.clone(),
1235 sources: polygon.sources.clone(),
1236 }];
1237 let mut convex: Vec<Cell> = Vec::new();
1238 let mut rounds = 0usize;
1239 while let Some(cell) = pending.pop() {
1240 rounds += 1;
1241 if rounds > 256 {
1242 return Err("the reflex cut did not converge in 256 rounds".into());
1243 }
1244 let cell = dedupe(cell, band);
1245 if cell.points.len() < 3 || cell.signed_area().abs() <= area_floor {
1246 continue;
1247 }
1248 match cell.reflex_vertex(band) {
1249 None => convex.push(cell),
1250 Some(vertex) => {
1251 let (first, second) = split_at_reflex(&cell, vertex, band)?;
1252 pending.push(first);
1253 pending.push(second);
1254 }
1255 }
1256 }
1257 let Some(hinge) = hinge else { return Ok(convex) };
1258 let mut cells: Vec<Cell> = Vec::new();
1259 for cell in convex {
1260 let sides: Vec<f64> = cell.points.iter().map(|point| hinge.side(*point)).collect();
1261 let high = sides.iter().copied().fold(f64::MIN, f64::max);
1262 let low = sides.iter().copied().fold(f64::MAX, f64::min);
1263 if high > band && low < -band {
1264 for half in split_by_line(&cell, &sides, band) {
1265 let half = dedupe(half, band);
1266 if half.points.len() >= 3 && half.signed_area().abs() > area_floor {
1267 cells.push(half);
1268 }
1269 }
1270 } else {
1271 cells.push(cell);
1272 }
1273 }
1274 Ok(cells)
1275}
1276
1277fn dedupe(cell: Cell, band: f64) -> Cell {
1280 let mut points: Vec<[f64; 2]> = Vec::new();
1281 let mut sources: Vec<EdgeSource> = Vec::new();
1282 for (point, source) in cell.points.iter().zip(&cell.sources) {
1283 if let Some(last) = points.last() {
1284 if (last[0] - point[0]).hypot(last[1] - point[1]) <= band {
1285 *sources.last_mut().expect("as many sources as points") = *source;
1287 continue;
1288 }
1289 }
1290 points.push(*point);
1291 sources.push(*source);
1292 }
1293 while points.len() >= 2 {
1294 let first = points[0];
1295 let last = points[points.len() - 1];
1296 if (last[0] - first[0]).hypot(last[1] - first[1]) <= band {
1297 points.pop();
1298 let source = sources.pop().expect("as many sources as points");
1299 let _ = source;
1300 } else {
1301 break;
1302 }
1303 }
1304 Cell { points, sources }
1305}
1306
1307fn split_at_reflex(cell: &Cell, vertex: usize, band: f64) -> Result<(Cell, Cell), String> {
1312 let n = cell.points.len();
1313 let points: Vec<[f64; 2]> = (0..n).map(|i| cell.points[(vertex + i) % n]).collect();
1315 let sources: Vec<EdgeSource> = (0..n).map(|i| cell.sources[(vertex + i) % n]).collect();
1316 let p0 = points[0];
1317 let prev = points[n - 1];
1318 let d = [p0[0] - prev[0], p0[1] - prev[1]];
1319 let length = (d[0] * d[0] + d[1] * d[1]).sqrt();
1320 if length <= band {
1321 return Err("a reflex corner with a zero-length arriving edge".into());
1322 }
1323 let d = [d[0] / length, d[1] / length];
1324 let mut best: Option<(usize, f64, f64)> = None;
1326 for j in 1..(n - 1) {
1327 let a = points[j];
1328 let b = points[(j + 1) % n];
1329 let e = [b[0] - a[0], b[1] - a[1]];
1330 let denominator = d[0] * e[1] - d[1] * e[0];
1331 if denominator.abs() <= 1e-15 {
1332 continue;
1333 }
1334 let w = [a[0] - p0[0], a[1] - p0[1]];
1335 let t = (w[0] * e[1] - w[1] * e[0]) / denominator;
1336 let s = (w[0] * d[1] - w[1] * d[0]) / denominator;
1337 let e_length = (e[0] * e[0] + e[1] * e[1]).sqrt().max(1e-300);
1338 let s_band = band / e_length;
1339 if t <= band || s < -s_band || s > 1.0 + s_band {
1340 continue;
1341 }
1342 if best.map_or(true, |(_, best_t, _)| t < best_t) {
1343 best = Some((j, t, s.clamp(0.0, 1.0)));
1344 }
1345 }
1346 let Some((j, t, s)) = best else {
1347 return Err(format!("the extension of the edge arriving at corner {vertex} meets no boundary edge"));
1348 };
1349 let hit = [p0[0] + d[0] * t, p0[1] + d[1] * t];
1350 let a = points[j];
1351 let b = points[(j + 1) % n];
1352 let e_length = ((b[0] - a[0]).powi(2) + (b[1] - a[1]).powi(2)).sqrt().max(1e-300);
1353 let at_a = s * e_length <= band;
1354 let at_b = (1.0 - s) * e_length <= band;
1355 let arriving = sources[n - 1];
1356 let mut first_points: Vec<[f64; 2]> = points[0..=j].to_vec();
1358 let mut first_sources: Vec<EdgeSource> = sources[0..j].to_vec();
1359 if at_b {
1360 first_points.push(b);
1361 first_sources.push(sources[j]);
1362 } else if !at_a {
1363 first_points.push(hit);
1364 first_sources.push(sources[j]);
1365 }
1366 first_sources.push(arriving);
1367 let mut second_points: Vec<[f64; 2]> = Vec::new();
1370 let mut second_sources: Vec<EdgeSource> = Vec::new();
1371 if at_a {
1372 second_points.push(a);
1373 second_sources.push(sources[j]);
1374 } else if !at_b {
1375 second_points.push(hit);
1376 second_sources.push(sources[j]);
1377 }
1378 for k in (j + 1)..n {
1379 second_points.push(points[k]);
1380 second_sources.push(sources[k]);
1381 }
1382 *second_sources.last_mut().expect("the second piece has an edge") = arriving;
1384 Ok((
1385 Cell { points: first_points, sources: first_sources },
1386 Cell { points: second_points, sources: second_sources },
1387 ))
1388}
1389
1390fn split_by_line(cell: &Cell, sides: &[f64], band: f64) -> [Cell; 2] {
1393 let n = cell.points.len();
1394 let mut positive = Cell { points: Vec::new(), sources: Vec::new() };
1395 let mut negative = Cell { points: Vec::new(), sources: Vec::new() };
1396 for j in 0..n {
1397 let p = cell.points[j];
1398 let q = cell.points[(j + 1) % n];
1399 let (fp, fq) = (sides[j], sides[(j + 1) % n]);
1400 let source = cell.sources[j];
1401 let p_on = fp.abs() <= band;
1402 let q_on = fq.abs() <= band;
1403 if fp > band || p_on {
1408 positive.points.push(p);
1409 positive.sources.push(if !p_on || fq > band || q_on { source } else { EdgeSource::Hinge });
1410 }
1411 if fp < -band || p_on {
1412 negative.points.push(p);
1413 negative.sources.push(if !p_on || fq < -band || q_on { source } else { EdgeSource::Hinge });
1414 }
1415 if (fp > band && fq < -band) || (fp < -band && fq > band) {
1417 let t = fp / (fp - fq);
1418 let x = [p[0] + (q[0] - p[0]) * t, p[1] + (q[1] - p[1]) * t];
1419 if fp > band {
1420 positive.points.push(x);
1421 positive.sources.push(EdgeSource::Hinge);
1422 negative.points.push(x);
1423 negative.sources.push(source);
1424 } else {
1425 negative.points.push(x);
1426 negative.sources.push(EdgeSource::Hinge);
1427 positive.points.push(x);
1428 positive.sources.push(source);
1429 }
1430 }
1431 }
1432 [positive, negative]
1433}
1434
1435fn one_way_only(
1443 id: u64,
1444 clip: &[HalfSpace],
1445 sides: &[HalfSpace],
1446 other: [HalfSpace; 2],
1447 tolerance: f64,
1448 name: impl Fn() -> String,
1449) -> Result<(), KernelRefusal> {
1450 let band = tolerance.max(1e-9) * 10.0;
1451 let spaces: Vec<HalfSpace> = clip.iter().chain(sides).chain(&other).copied().collect();
1452 let mut depth = 0.0f64;
1453 for a in 0..spaces.len() {
1454 for b in (a + 1)..spaces.len() {
1455 for c in (b + 1)..spaces.len() {
1456 let Some(point) = three_plane_point(&spaces[a], &spaces[b], &spaces[c]) else {
1457 continue;
1458 };
1459 if spaces
1460 .iter()
1461 .any(|space| space.normal.dot(point) - space.offset > band)
1462 {
1463 continue;
1464 }
1465 let inside = other
1467 .iter()
1468 .map(|space| space.offset - space.normal.dot(point))
1469 .fold(0.0f64, f64::max);
1470 depth = depth.max(inside);
1471 }
1472 }
1473 }
1474 if depth > band * 1e3 {
1475 return Err(KernelRefusal::ill_posed(
1476 KernelStage::Classify,
1477 "one_way_only",
1478 format!(
1479 "{OP}: the motion turns face {id}{} so that the line its old and new carriers meet \
1480 in crosses its own bounds — part of the face moves into the body and part out of it \
1481 (the other wedge reaches {depth:.3e} past them), so one tool would add material \
1482 where the other takes it away; refusing rather than choosing",
1483 name()
1484 ),
1485 ));
1486 }
1487 Ok(())
1488}
1489
1490fn cell_is_its_bounds(
1505 face: &FaceRecord,
1506 cell: &Cell,
1507 frame: &Frame,
1508 at_home: &[HalfSpace],
1509 scale: f64,
1510) -> Result<(), KernelRefusal> {
1511 let band = (scale * 1e-4).max(1e-6);
1516 let n = cell.points.len();
1517 let mut worst = 0.0f64;
1518 for i in 0..n {
1519 let a = cell.points[i];
1520 let b = cell.points[(i + 1) % n];
1521 for point in [a, [(a[0] + b[0]) * 0.5, (a[1] + b[1]) * 0.5]] {
1522 let point = frame.lift(point);
1523 for space in at_home {
1524 worst = worst.max(space.normal.dot(point) - space.offset);
1525 }
1526 }
1527 }
1528 if worst > band {
1529 return Err(KernelRefusal::unsupported(
1530 KernelStage::Classify,
1531 "reflex_trim",
1532 format!(
1533 "{OP}: face {}{} has a boundary point {worst:.3e} outside one of its own \
1534 bounding carriers, so its trim is not the region its bounding carriers \
1535 enclose — their intersection is the trim's CONVEX CORE, and the tool this \
1536 road would build is that core rather than all the material the motion \
1537 moves. Refusing rather than cutting some of it",
1538 face.id,
1539 named(face)
1540 ),
1541 ));
1542 }
1543 Ok(())
1544}
1545
1546
1547fn face_by_id(solid: &BrepSolid, id: u64) -> Result<&FaceRecord, KernelRefusal> {
1548 solid
1549 .shells
1550 .iter()
1551 .flat_map(|shell| &shell.faces)
1552 .find(|face| face.id == id)
1553 .ok_or_else(|| {
1554 KernelRefusal::input(
1555 KernelStage::Collect,
1556 "face_id",
1557 format!("{OP}: no face with id {id} on this solid"),
1558 )
1559 })
1560}
1561
1562fn named(face: &FaceRecord) -> String {
1563 face.name
1564 .as_ref()
1565 .map(|name| format!(" '{name}'"))
1566 .unwrap_or_default()
1567}
1568
1569fn solid_bounds(solid: &BrepSolid) -> Result<(Vec3, Vec3), KernelRefusal> {
1573 if solid.vertices.is_empty() {
1574 return Err(KernelRefusal::input(
1575 KernelStage::Collect,
1576 "empty_solid",
1577 format!("{OP}: the solid has no vertices to bound"),
1578 ));
1579 }
1580 let mut low = Vec3::new(f64::MAX, f64::MAX, f64::MAX);
1581 let mut high = Vec3::new(f64::MIN, f64::MIN, f64::MIN);
1582 for vertex in &solid.vertices {
1583 low = Vec3::new(
1584 low.x.min(vertex.point.x),
1585 low.y.min(vertex.point.y),
1586 low.z.min(vertex.point.z),
1587 );
1588 high = Vec3::new(
1589 high.x.max(vertex.point.x),
1590 high.y.max(vertex.point.y),
1591 high.z.max(vertex.point.z),
1592 );
1593 }
1594 Ok((low, high))
1595}
1596
1597fn box_half_spaces((low, high): (Vec3, Vec3), margin: f64) -> Vec<HalfSpace> {
1599 let low = Vec3::new(low.x - margin, low.y - margin, low.z - margin);
1600 let high = Vec3::new(high.x + margin, high.y + margin, high.z + margin);
1601 vec![
1602 HalfSpace { normal: Vec3::new(1.0, 0.0, 0.0), offset: high.x },
1603 HalfSpace { normal: Vec3::new(-1.0, 0.0, 0.0), offset: -low.x },
1604 HalfSpace { normal: Vec3::new(0.0, 1.0, 0.0), offset: high.y },
1605 HalfSpace { normal: Vec3::new(0.0, -1.0, 0.0), offset: -low.y },
1606 HalfSpace { normal: Vec3::new(0.0, 0.0, 1.0), offset: high.z },
1607 HalfSpace { normal: Vec3::new(0.0, 0.0, -1.0), offset: -low.z },
1608 ]
1609}
1610
1611fn convex_solid(
1622 spaces: &[HalfSpace],
1623 inside: Vec3,
1624 tolerance: f64,
1625) -> Result<BrepSolid, String> {
1626 let band = tolerance.max(1e-9) * 10.0;
1627 for space in spaces {
1628 if space.normal.dot(inside) - space.offset > -band {
1629 return Err(format!(
1630 "the seat point is not strictly inside the half-space set (it is {:.3e} past one \
1631 of them), so the slab and its bounds enclose nothing",
1632 space.normal.dot(inside) - space.offset
1633 ));
1634 }
1635 }
1636 let mut corners: Vec<Vec3> = Vec::new();
1638 for a in 0..spaces.len() {
1639 for b in (a + 1)..spaces.len() {
1640 for c in (b + 1)..spaces.len() {
1641 let Some(point) = three_plane_point(&spaces[a], &spaces[b], &spaces[c]) else {
1642 continue;
1643 };
1644 if spaces
1645 .iter()
1646 .any(|space| space.normal.dot(point) - space.offset > band)
1647 {
1648 continue;
1649 }
1650 if corners
1651 .iter()
1652 .any(|other| other.sub(point).length() <= band)
1653 {
1654 continue;
1655 }
1656 corners.push(point);
1657 }
1658 }
1659 }
1660 if corners.len() < 4 {
1661 return Err(format!(
1662 "the half-space set has {} corner(s), which is not a bounded volume",
1663 corners.len()
1664 ));
1665 }
1666
1667 let vertices: Vec<VertexRecord> = corners
1668 .iter()
1669 .enumerate()
1670 .map(|(index, point)| VertexRecord { id: index as u64 + 1, point: *point })
1671 .collect();
1672 let mut edges: Vec<EdgeRecord> = Vec::new();
1673 let mut edge_index: HashMap<(usize, usize), u64> = HashMap::default();
1674 let mut faces: Vec<FaceRecord> = Vec::new();
1675 let mut next_id = 1000u64;
1676
1677 for space in spaces {
1678 let mut on: Vec<usize> = (0..corners.len())
1679 .filter(|&index| (space.normal.dot(corners[index]) - space.offset).abs() <= band)
1680 .collect();
1681 if on.len() < 3 {
1682 continue;
1683 }
1684 let centroid = on
1688 .iter()
1689 .fold(Vec3::new(0.0, 0.0, 0.0), |sum, &index| sum.add(corners[index]))
1690 .scale(1.0 / on.len() as f64);
1691 let u_dir = corners[on[0]].sub(centroid).normalized()?;
1692 let v_dir = space.normal.cross(u_dir).normalized()?;
1693 on.sort_by(|&a, &b| {
1694 let angle = |index: usize| {
1695 let delta = corners[index].sub(centroid);
1696 delta.dot(v_dir).atan2(delta.dot(u_dir))
1697 };
1698 angle(a).total_cmp(&angle(b))
1699 });
1700 faces.push(polygon_face(
1701 &vertices,
1702 &on,
1703 space.normal,
1704 &mut edges,
1705 &mut edge_index,
1706 &mut next_id,
1707 )?);
1708 }
1709 if faces.len() < 4 {
1710 return Err(format!(
1711 "the half-space set produced {} face(s), which does not close a volume",
1712 faces.len()
1713 ));
1714 }
1715
1716 let shell_id = next_id + 1;
1717 let solid_id = next_id + 2;
1718 let solid = BrepSolid {
1719 mass_properties_cache: Default::default(),
1720 cells: Default::default(),
1721 id: solid_id,
1722 vertices,
1723 edges,
1724 shells: vec![ShellRecord { id: shell_id, faces }],
1725 genus: 0,
1726 kind: crate::BodyKind::Solid,
1727 };
1728 let issues = solid.validate();
1729 if !issues.is_empty() {
1730 return Err(format!("the tool polytope did not validate: {issues:?}"));
1731 }
1732 Ok(solid)
1733}
1734
1735fn three_plane_point(a: &HalfSpace, b: &HalfSpace, c: &HalfSpace) -> Option<Vec3> {
1738 let cross = b.normal.cross(c.normal);
1739 let determinant = a.normal.dot(cross);
1740 if determinant.abs() <= 1e-9 {
1741 return None;
1742 }
1743 let point = cross
1744 .scale(a.offset)
1745 .add(c.normal.cross(a.normal).scale(b.offset))
1746 .add(a.normal.cross(b.normal).scale(c.offset))
1747 .scale(1.0 / determinant);
1748 Some(point)
1749}
1750
1751fn polygon_face(
1754 vertices: &[VertexRecord],
1755 corners: &[usize],
1756 normal: Vec3,
1757 edges: &mut Vec<EdgeRecord>,
1758 edge_index: &mut HashMap<(usize, usize), u64>,
1759 next_id: &mut u64,
1760) -> Result<FaceRecord, String> {
1761 let origin_point = vertices[corners[0]].point;
1762 let u_dir = vertices[corners[1]]
1763 .point
1764 .sub(origin_point)
1765 .normalized()?;
1766 let v_dir = normal.cross(u_dir).normalized()?;
1767 let mut low_u = f64::MAX;
1768 let mut low_v = f64::MAX;
1769 let mut high_u = f64::MIN;
1770 let mut high_v = f64::MIN;
1771 let projected: Vec<(f64, f64)> = corners
1772 .iter()
1773 .map(|&index| {
1774 let delta = vertices[index].point.sub(origin_point);
1775 let (u, v) = (delta.dot(u_dir), delta.dot(v_dir));
1776 low_u = low_u.min(u);
1777 low_v = low_v.min(v);
1778 high_u = high_u.max(u);
1779 high_v = high_v.max(v);
1780 (u, v)
1781 })
1782 .collect();
1783 let plane_origin = origin_point.add(u_dir.scale(low_u)).add(v_dir.scale(low_v));
1784 let surface = make_plane(
1785 plane_origin,
1786 u_dir,
1787 v_dir,
1788 high_u - low_u,
1789 high_v - low_v,
1790 )?;
1791
1792 let mut coedges: Vec<CoedgeRecord> = Vec::with_capacity(corners.len());
1793 for index in 0..corners.len() {
1794 let start = corners[index];
1795 let end = corners[(index + 1) % corners.len()];
1796 let key = if start < end { (start, end) } else { (end, start) };
1797 let edge_id = match edge_index.get(&key) {
1798 Some(id) => *id,
1799 None => {
1800 let id = *next_id;
1801 *next_id += 1;
1802 edges.push(EdgeRecord {
1803 id,
1804 curve: make_line(vertices[key.0].point, vertices[key.1].point)?,
1805 t0: 0.0,
1806 t1: 1.0,
1807 start_vertex_id: vertices[key.0].id,
1808 end_vertex_id: vertices[key.1].id,
1809 degenerate: false,
1810 name: None,
1811 });
1812 edge_index.insert(key, id);
1813 id
1814 }
1815 };
1816 let start_uv = projected[index];
1817 let end_uv = projected[(index + 1) % corners.len()];
1818 coedges.push(CoedgeRecord {
1819 id: *next_id,
1820 edge_id,
1821 forward: vertices[start].id
1822 == edges
1823 .iter()
1824 .find(|edge| edge.id == edge_id)
1825 .expect("the edge was just minted")
1826 .start_vertex_id,
1827 pcurve: parameter_line(
1828 start_uv.0 - low_u,
1829 start_uv.1 - low_v,
1830 end_uv.0 - low_u,
1831 end_uv.1 - low_v,
1832 )?,
1833 });
1834 *next_id += 1;
1835 }
1836 let loop_id = *next_id;
1837 *next_id += 1;
1838 let face_id = *next_id;
1839 *next_id += 1;
1840 Ok(FaceRecord {
1841 id: face_id,
1842 surface,
1843 same_sense: true,
1844 loops: vec![LoopRecord { id: loop_id, coedges }],
1845 name: None,
1846 })
1847}
1848
1849fn restamp(
1859 mut cut: BrepSolid,
1860 original: &BrepSolid,
1861 moved: &HashSet<u64>,
1862 motion: &Motion,
1863 tolerance: f64,
1864) -> BrepSolid {
1865 let band = tolerance.max(1e-9) * 1e3;
1866 let mut wanted: Vec<(String, Vec3, f64)> = Vec::new();
1868 for face in original.shells.iter().flat_map(|shell| &shell.faces) {
1869 let Some(name) = face.name.clone() else { continue };
1870 let Ok(plane) = plane_of_surface(&face.surface, band, OP) else {
1871 continue;
1872 };
1873 let normal = if face.same_sense {
1874 plane.normal
1875 } else {
1876 plane.normal.scale(-1.0)
1877 };
1878 let (normal, offset) = if moved.contains(&face.id) {
1879 motion.plane(normal, normal.dot(plane.origin))
1880 } else {
1881 (normal, normal.dot(plane.origin))
1882 };
1883 wanted.push((name, normal, offset));
1884 }
1885 let taken: HashSet<String> = cut
1886 .shells
1887 .iter()
1888 .flat_map(|shell| &shell.faces)
1889 .filter_map(|face| face.name.clone())
1890 .collect();
1891 let mut used: HashMap<String, usize> = HashMap::default();
1892 for name in &taken {
1893 used.insert(name.clone(), 1);
1894 }
1895 for shell in &mut cut.shells {
1896 for face in &mut shell.faces {
1897 if face.name.is_some() {
1898 continue;
1899 }
1900 let Ok(plane) = plane_of_surface(&face.surface, band, OP) else {
1901 continue;
1902 };
1903 let normal = if face.same_sense {
1904 plane.normal
1905 } else {
1906 plane.normal.scale(-1.0)
1907 };
1908 let offset = normal.dot(plane.origin);
1909 let Some((name, _, _)) = wanted.iter().find(|(_, other, other_offset)| {
1910 other.sub(normal).length() <= 1e-9 && (other_offset - offset).abs() <= band
1911 }) else {
1912 continue;
1913 };
1914 let count = used.entry(name.clone()).or_insert(0);
1915 *count += 1;
1916 face.name = Some(if *count == 1 {
1917 name.clone()
1918 } else {
1919 format!("{name}_{}", *count - 1)
1920 });
1921 }
1922 }
1923 cut
1924}