geop_lib_booleans/
split.rs

1//! Splitting a solid into pieces with a sheet.
2//!
3//! The sheet is imprinted onto the solid exactly as a boolean imprints two
4//! solids ([`remesh`]), after which every piece of the sheet lies wholly
5//! inside the solid or wholly outside it. The pieces inside are where the
6//! solid is cut. Each becomes a face of *two* pieces of the solid: as it is
7//! on the side its normal points away from, turned around on the other.
8//!
9//! Which faces bound one piece is then purely topological. A cut face meets
10//! the solid's boundary along the edges where the sheet crossed it, with one
11//! face of the solid on either side; of the two, the one running along the
12//! edge the other way than the cut face does is the one on the cut face's
13//! side — two faces of a closed, oriented shell always run their shared edge
14//! in opposite directions. Elsewhere, faces sharing an edge bound the same
15//! piece. So the pieces are the connected sets of "face sides", and each is
16//! built as a solid of its own: the cut sides it shares with its neighbours
17//! are copies, so the pieces share nothing.
18
19use std::collections::{HashMap, HashSet};
20
21use geop_core_math::{
22    geop_error::{GeopError, GeopResult, WithContext},
23    union_find::UnionFind,
24    with_context,
25};
26use geop_core_topology::{
27    Body, CoedgeGeometry, EdgeId, FaceId, ShellId, SolidId,
28    build::{BodySpec, BuiltBody},
29    contains::shell::Classified,
30};
31use geop_ops_model::{BodyNames, NamedModel, Namer, RefId};
32
33use crate::{
34    boolean::{FaceClassification, classify_face},
35    remesh::{RemeshParams, remesh},
36};
37
38/// Splits `solid` along the sheet `sheet` into the pieces it cuts the solid
39/// into, each a solid of its own. The solid is consumed; the sheet is left as
40/// it is — a copy of it does the cutting.
41///
42/// The pieces are named `N(0)`, `N(1)`, ..., ordered by the smallest name of
43/// a face they kept of the solid, so the order does not depend on how they
44/// were found. Every face, edge and vertex keeps its name in the first piece
45/// it is part of; a piece sharing it with an earlier one — the cut, on its
46/// other side, and the edges and vertices around it — names its own copy
47/// `N(name,k)` for piece `k`. What the cutting creates is named as a
48/// boolean names what it creates (see `naming`), the copy of the
49/// sheet's face `F` being `N(F)`.
50///
51/// Fails, leaving the solid in pieces of no use, unless the sheet cuts the
52/// solid clean through: a sheet ending inside the solid, or one along which
53/// the solid stays connected, splits nothing.
54pub fn split(
55    model: &mut NamedModel,
56    namer: &Namer,
57    solid: SolidId,
58    sheet: ShellId,
59    params: RemeshParams,
60) -> GeopResult<Vec<SolidId>> {
61    let ctx = with_context!("split(name={}, solid={solid}, sheet={sheet})", namer.root());
62
63    // The cutting copy of the sheet.
64    let sheet_faces = model
65        .topology()
66        .body_faces(Body::Sheet(sheet))
67        .with_context(&ctx)?;
68    let cutter = model
69        .copy_faces(&sheet_faces, None, |name| namer.name(&[name]))
70        .and_then(|built| built.body())
71        .with_context(&ctx)?;
72
73    remesh(model, namer, solid, cutter, params).with_context(&ctx)?;
74    let topology = model.topology();
75    let solid_faces = topology.solid_faces(solid).with_context(&ctx)?;
76    let mut cuts = Vec::new();
77    let mut classified = Classified::default();
78    for face in topology.body_faces(cutter).with_context(&ctx)? {
79        let class = classify_face(topology, face, solid, params, &mut classified)
80            .with_context(&ctx)
81            .with_context(with_context!("classifying face {face}"))?;
82        if class == FaceClassification::Inside {
83            cuts.push(face);
84        }
85    }
86    if cuts.is_empty() {
87        return Err(ctx(GeopError::new(
88            "split: the face does not cut through the solid",
89        )));
90    }
91
92    // Face sides: each face of the solid, then each cut as it is ("back",
93    // the side its normal points away from), then each turned around.
94    let n = solid_faces.len();
95    let k = cuts.len();
96    let side: HashMap<FaceId, usize> = solid_faces
97        .iter()
98        .chain(&cuts)
99        .enumerate()
100        .map(|(i, &f)| (f, i))
101        .collect();
102    let turned = |i: usize| i + k;
103    let mut pieces = UnionFind::new(n + 2 * k);
104    let mut users: HashMap<EdgeId, Vec<(usize, geop_core_topology::Sense)>> = HashMap::new();
105    // In the order of the sides, not of the map: it decides the order of
106    // each edge's users below.
107    for (i, &face) in solid_faces.iter().chain(&cuts).enumerate() {
108        for coedge in topology.iterate_face_coedges(face) {
109            let coedge = topology.get_coedge(coedge)?;
110            if let CoedgeGeometry::Edge(edge) = coedge.geometry {
111                users.entry(edge).or_default().push((i, coedge.sense));
112            }
113        }
114    }
115    let sides_of = |i: usize| if i < n { solid_faces[i] } else { cuts[i - n] };
116    let mut edges: Vec<_> = users.into_iter().collect();
117    edges.sort_by_key(|(e, _)| e.0);
118    for (edge, uses) in edges {
119        let (on_solid, on_cut): (Vec<_>, Vec<_>) = uses.into_iter().partition(|&(i, _)| i < n);
120        match (&on_solid[..], &on_cut[..]) {
121            ([(a, _), (b, _)], []) => pieces.union(*a, *b),
122            // The sheet crossed the solid's boundary here.
123            ([(a, sense_a), (b, _)], [(cut, sense)]) => {
124                let (with, against) = if sense_a != sense { (*a, *b) } else { (*b, *a) };
125                pieces.union(*cut, with);
126                pieces.union(turned(*cut), against);
127            }
128            // Two cuts meeting inside the solid.
129            ([], [(c, sense_c), (d, sense_d)]) => {
130                if sense_c != sense_d {
131                    pieces.union(*c, *d);
132                    pieces.union(turned(*c), turned(*d));
133                } else {
134                    pieces.union(*c, turned(*d));
135                    pieces.union(turned(*c), *d);
136                }
137            }
138            ([], [(cut, _)]) => {
139                let users: Vec<_> = (topology.coedges_of_edge(edge).iter())
140                    .filter_map(|&c| topology.get_coedge(c).ok())
141                    .map(|c| model.name_of(c.face).unwrap_or("?").to_string())
142                    .collect();
143                return Err(ctx(GeopError::new(format!(
144                    "split: the face ends inside the solid, at edge {edge} ({:?}) of the cut {:?}, used by the faces {users:?}; it must cut all the way through",
145                    model.name_of(edge),
146                    model.name_of(sides_of(*cut)),
147                ))));
148            }
149            _ => {
150                return Err(ctx(GeopError::new(format!(
151                    "split: {} face(s) of the solid and {} of the cut meet at edge {edge}, which does not divide the solid into pieces",
152                    on_solid.len(),
153                    on_cut.len()
154                ))));
155            }
156        }
157    }
158
159    let groups = pieces.groups();
160    let piece_of: HashMap<usize, usize> = groups
161        .iter()
162        .enumerate()
163        .flat_map(|(g, members)| members.iter().map(move |&m| (m, g)))
164        .collect();
165    if let Some(&cut) = cuts
166        .iter()
167        .find(|c| piece_of[&side[c]] == piece_of[&turned(side[c])])
168    {
169        return Err(ctx(GeopError::new(format!(
170            "split: the solid stays connected around the cut face {cut}: the face does not divide it into pieces"
171        ))));
172    }
173
174    // Each piece, described: its faces of the solid as they are, its cuts as
175    // they are or turned around.
176    let sides: Vec<FaceId> = solid_faces
177        .iter()
178        .chain(&cuts)
179        .chain(&cuts)
180        .copied()
181        .collect();
182    let mut described: Vec<(String, BodySpec, Vec<RefId>)> = Vec::new();
183    for members in &groups {
184        let faces: Vec<FaceId> = members.iter().map(|&m| sides[m]).collect();
185        let (mut spec, sources) = topology.body_spec(&faces, true).with_context(&ctx)?;
186        for (face, &m) in spec.faces.iter_mut().zip(members) {
187            if m >= n + k {
188                *face = face.reversed();
189            }
190        }
191        let first = members
192            .iter()
193            .filter(|&&m| m < n)
194            .map(|&m| model.name(solid_faces[m]).map(str::to_string))
195            .collect::<GeopResult<Vec<_>>>()?
196            .into_iter()
197            .min()
198            .ok_or_else(|| ctx(GeopError::new("split: a piece has no face of the solid")))?;
199        let entities = sources
200            .vertices
201            .iter()
202            .map(|&v| RefId::from(v))
203            .chain(sources.edges.iter().map(|&e| e.into()))
204            .chain(sources.faces.iter().map(|&f| f.into()))
205            .collect();
206        described.push((first, spec, entities));
207    }
208    described.sort_by(|a, b| a.0.cmp(&b.0));
209
210    // Names, before the originals and their names are gone.
211    let mut taken: HashSet<RefId> = HashSet::new();
212    let mut named: Vec<(BodySpec, BodyNames)> = Vec::new();
213    for (index, (_, spec, entities)) in described.into_iter().enumerate() {
214        let mut all = Vec::with_capacity(entities.len());
215        for id in &entities {
216            let original = model.name(*id).map(str::to_string)?;
217            all.push(if taken.insert(*id) {
218                original
219            } else {
220                namer.name(&[&original, &index.to_string()])
221            });
222        }
223        let (v, e) = (spec.vertices.len(), spec.edges.len());
224        let names = BodyNames {
225            vertices: all[..v].to_vec(),
226            edges: all[v..v + e].to_vec(),
227            faces: all[v + e..].to_vec(),
228            solid: Some(namer.name(&[&index.to_string()])),
229        };
230        named.push((spec, names));
231    }
232
233    model
234        .assemble_solid(&[solid.into(), cutter], &[], namer.root())
235        .with_context(&ctx)?;
236    named
237        .into_iter()
238        .map(|(spec, names)| {
239            let BuiltBody { solid, .. } = model.build_body(spec, names).with_context(&ctx)?;
240            Ok(solid.expect("built as a solid"))
241        })
242        .collect()
243}
244
245#[cfg(test)]
246mod tests {
247    use super::split;
248    use crate::remesh::RemeshParams;
249    use geop_core_math::{primitives::CoordinateSystem, scalars::Scalar, vector::Vector3};
250    use geop_core_topology::{
251        Body,
252        contains::shell::{Classified, PointClassification, solid_contains},
253    };
254    use geop_lib_sweep::{
255        SweepLoop, extrude, polyline,
256        shapes::{block, v2, v3},
257    };
258    use geop_ops_model::Profile;
259    use geop_ops_model::{NamedModel, Namer, testing::assert_valid};
260
261    /// A unit cube around the origin, and a sheet: the chain through
262    /// `points` in the `z = -1` plane, swept up to `z = 1`.
263    fn cube_and_sheet(points: &[[f64; 2]]) -> (NamedModel, geop_core_topology::SolidId, Body) {
264        let mut model = NamedModel::new();
265        let cube = block(&mut model, "c", [-0.5; 3], [0.5; 3]);
266        let chain = polyline(&points.iter().map(|p| v2(p[0], p[1])).collect::<Vec<_>>()).unwrap();
267        let plane = CoordinateSystem::world_at(v3(0.0, 0.0, -1.0));
268        let namer = Namer::new("extrude", "cut").unwrap();
269        let built = extrude(
270            &mut model,
271            &namer,
272            None,
273            &plane,
274            &[Scalar::ZERO, Scalar::from_f64(2.0)],
275            &[SweepLoop::plain(Profile::open(chain))],
276        )
277        .unwrap();
278        (model, cube, built.body().unwrap())
279    }
280
281    /// Which of `pieces` holds `p`.
282    fn holding(
283        model: &NamedModel,
284        pieces: &[geop_core_topology::SolidId],
285        p: Vector3,
286    ) -> Vec<usize> {
287        (0..pieces.len())
288            .filter(|&k| {
289                solid_contains(
290                    model.topology(),
291                    pieces[k],
292                    p,
293                    2000,
294                    Scalar::from_f64(1e-7),
295                    7,
296                    &mut Classified::default(),
297                )
298                .unwrap()
299                    == PointClassification::Inside
300            })
301            .collect()
302    }
303
304    /// A flat sheet across the cube cuts it in two halves, and stays as it
305    /// was.
306    #[test]
307    fn a_flat_sheet_cuts_a_cube_in_two() {
308        let (mut model, cube, sheet) = cube_and_sheet(&[[0.1, -1.0], [0.1, 1.0]]);
309        let Body::Sheet(shell) = sheet else {
310            unreachable!()
311        };
312        let namer = Namer::new("split", "s").unwrap();
313        let pieces = split(&mut model, &namer, cube, shell, RemeshParams::default()).unwrap();
314        assert_valid(&model);
315        assert_eq!(pieces.len(), 2);
316        assert_eq!(model.solid_names(), ["split(s,0)", "split(s,1)"]);
317        let (left, right) = (
318            holding(&model, &pieces, v3(-0.2, 0.0, 0.0)),
319            holding(&model, &pieces, v3(0.3, 0.0, 0.0)),
320        );
321        assert_eq!(left.len(), 1);
322        assert_eq!(right.len(), 1);
323        assert_ne!(left, right);
324        // Each half has the cube's six faces, cut down, the right one five
325        // of them and its own side of the cut.
326        let faces = |k: usize| model.topology().solid_faces(pieces[k]).unwrap().len();
327        assert_eq!(faces(0) + faces(1), 6 + 4 + 2);
328        assert_eq!(model.topology().body_faces(sheet).unwrap().len(), 1);
329    }
330
331    /// A sheet bent into a "U" whose bottom lies in the cube cuts the slot
332    /// inside the "U" out of it; the rest stays one piece, around the
333    /// bottom.
334    #[test]
335    fn a_bent_sheet_cuts_a_slot_out_of_a_cube() {
336        let (mut model, cube, sheet) =
337            cube_and_sheet(&[[-0.2, 1.0], [-0.2, -0.2], [0.2, -0.2], [0.2, 1.0]]);
338        let Body::Sheet(shell) = sheet else {
339            unreachable!()
340        };
341        let namer = Namer::new("split", "s").unwrap();
342        let pieces = split(&mut model, &namer, cube, shell, RemeshParams::default()).unwrap();
343        assert_valid(&model);
344        assert_eq!(pieces.len(), 2);
345        let slot = holding(&model, &pieces, v3(0.0, 0.2, 0.0));
346        assert_eq!(slot.len(), 1);
347        for p in [v3(-0.4, 0.0, 0.0), v3(0.4, 0.0, 0.0), v3(0.0, -0.4, 0.0)] {
348            let rest = holding(&model, &pieces, p);
349            assert!(rest.len() == 1 && rest != slot, "{p:?}");
350        }
351    }
352
353    /// A sheet missing the cube cuts nothing.
354    #[test]
355    fn a_sheet_beside_the_solid_cuts_nothing() {
356        let (mut model, cube, sheet) = cube_and_sheet(&[[2.0, -1.0], [2.0, 1.0]]);
357        let Body::Sheet(shell) = sheet else {
358            unreachable!()
359        };
360        let namer = Namer::new("split", "s").unwrap();
361        let error = split(&mut model, &namer, cube, shell, RemeshParams::default()).unwrap_err();
362        assert!(error.root_message().contains("does not cut"), "{error:?}");
363    }
364}
geop_lib_booleans/split.rs

geop_lib_booleans/
split.rs
