Kernel: pole-free six-patch spheres; vendored truck-shapeops with robust meshing and nearest-parameter loop projection (sphere and box booleans work); least-squares sphere classification

Oracle: 25/29 (14 and 15 now match Fusion). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01AS5Tn3xieewUSg3uxrMVdL
2026-10-07 21:30:18 +00:00
1 parent 7ef3426cb3
commit 043a542e29
@@ -0,0 +1,182 @@
+use truck_meshalgo::prelude::*;
+
+#[derive(Clone, Debug)]
+pub enum Alternative<T, U> {
+    FirstType(T),
+    SecondType(U),
+}
+
+#[doc(hidden)]
+#[macro_export]
+macro_rules! impl_from {
+    ($firsttype: ty, $secondtype: ty) => {
+        impl From<$firsttype> for $crate::alternative::Alternative<$firsttype, $secondtype> {
+            #[inline(always)]
+            fn from(t: $firsttype) -> Self { $crate::alternative::Alternative::FirstType(t) }
+        }
+    };
+}
+
+impl<T, U> From<U> for Alternative<T, U> {
+    #[inline(always)]
+    fn from(u: U) -> Self { Alternative::SecondType(u) }
+}
+
+// test for impl_from
+impl_from!((), usize);
+
+macro_rules! derive_method {
+	($method: tt, $return_type: ty, $($var: ident : $paramtype: ty),*) => {
+		fn $method (&self, $($var: $paramtype),*) -> $return_type {
+			match &self {
+				Alternative::FirstType(got) => got.$method($($var),*),
+				Alternative::SecondType(got) => got.$method($($var),*),
+			}
+		}
+	};
+	($method: tt <$x: ident : $y: path>, $return_type: ty, $($var: ident : $paramtype: ty),*) => {
+		fn $method<$x : $y>(&self, $($var: $paramtype),*) -> $return_type {
+			match &self {
+				Alternative::FirstType(got) => got.$method($($var),*),
+				Alternative::SecondType(got) => got.$method($($var),*),
+			}
+		}
+	};
+}
+
+impl<C0, C1> ParametricCurve for Alternative<C0, C1>
+where
+    C0: ParametricCurve,
+    C1: ParametricCurve<Point = C0::Point, Vector = C0::Vector>,
+{
+    type Point = C0::Point;
+    type Vector = C0::Vector;
+    derive_method!(subs, C0::Point, t: f64);
+    derive_method!(der, C0::Vector, t: f64);
+    derive_method!(der2, C0::Vector, t: f64);
+    derive_method!(
+        parameter_range,
+        (std::ops::Bound<f64>, std::ops::Bound<f64>),
+    );
+}
+
+impl<C0, C1> BoundedCurve for Alternative<C0, C1>
+where
+    C0: BoundedCurve,
+    C1: BoundedCurve<Point = C0::Point, Vector = C0::Vector>,
+{
+}
+
+impl<S0, S1> ParametricSurface for Alternative<S0, S1>
+where
+    S0: ParametricSurface,
+    S1: ParametricSurface<Point = S0::Point, Vector = S0::Vector>,
+{
+    type Point = S0::Point;
+    type Vector = S0::Vector;
+    derive_method!(subs, S0::Point, u: f64, v: f64);
+    derive_method!(uder, S0::Vector, u: f64, v: f64);
+    derive_method!(vder, S0::Vector, u: f64, v: f64);
+    derive_method!(uuder, S0::Vector, u: f64, v: f64);
+    derive_method!(uvder, S0::Vector, u: f64, v: f64);
+    derive_method!(vvder, S0::Vector, u: f64, v: f64);
+}
+
+impl<S0, S1> ParametricSurface3D for Alternative<S0, S1>
+where
+    S0: ParametricSurface3D,
+    S1: ParametricSurface3D,
+{
+    derive_method!(normal, Vector3, u: f64, v: f64);
+}
+
+impl<C0, C1> Cut for Alternative<C0, C1>
+where
+    C0: Cut,
+    C1: Cut<Point = C0::Point, Vector = C0::Vector>,
+{
+    fn cut(&mut self, t: f64) -> Self {
+        match self {
+            Self::FirstType(curve) => Self::FirstType(curve.cut(t)),
+            Self::SecondType(curve) => Self::SecondType(curve.cut(t)),
+        }
+    }
+}
+
+impl<C0, C1> ParameterDivision1D for Alternative<C0, C1>
+where
+    C0: ParameterDivision1D,
+    C1: ParameterDivision1D<Point = C0::Point>,
+{
+    type Point = C0::Point;
+    derive_method!(
+        parameter_division,
+        (Vec<f64>, Vec<C0::Point>),
+        range: (f64, f64),
+        tol: f64
+    );
+}
+
+impl<S0, S1> ParameterDivision2D for Alternative<S0, S1>
+where
+    S0: ParameterDivision2D,
+    S1: ParameterDivision2D,
+{
+    derive_method!(
+        parameter_division,
+        (Vec<f64>, Vec<f64>),
+        range: ((f64, f64), (f64, f64)),
+        tol: f64
+    );
+}
+
+impl<D: SPDimension, T, U> SearchParameter<D> for Alternative<T, U>
+where
+    T: SearchParameter<D>,
+    U: SearchParameter<D, Point = T::Point>,
+{
+    type Point = T::Point;
+    derive_method!(
+        search_parameter<H: Into<D::Hint>>,
+        Option<D::Parameter>,
+        point: T::Point,
+        hint: H,
+        trials: usize
+    );
+}
+
+impl<D: SPDimension, T, U> SearchNearestParameter<D> for Alternative<T, U>
+where
+    T: SearchNearestParameter<D>,
+    U: SearchNearestParameter<D, Point = T::Point>,
+{
+    type Point = T::Point;
+    derive_method!(
+        search_nearest_parameter<H: Into<D::Hint>>,
+        Option<D::Parameter>,
+        point: T::Point,
+        hint: H,
+        trials: usize
+    );
+}
+
+impl<T, U> Invertible for Alternative<T, U>
+where
+    T: Invertible,
+    U: Invertible,
+{
+    #[inline(always)]
+    fn invert(&mut self) {
+        match self {
+            Self::FirstType(entity) => entity.invert(),
+            Self::SecondType(entity) => entity.invert(),
+        }
+    }
+    #[inline(always)]
+    fn inverse(&self) -> Self {
+        match self {
+            Self::FirstType(entity) => Self::FirstType(entity.inverse()),
+            Self::SecondType(entity) => Self::SecondType(entity.inverse()),
+        }
+    }
+}
@@ -0,0 +1,132 @@
+use rustc_hash::{FxHashMap as HashMap, FxHashSet as HashSet};
+use truck_geometry::prelude::*;
+use truck_meshalgo::rexport_polymesh::*;
+use truck_topology::compress::*;
+
+type Edge<C> = CompressedEdge<C>;
+type EdgeIndex = CompressedEdgeIndex;
+type Wire = Vec<EdgeIndex>;
+type Face<S> = CompressedFace<S>;
+type Shell<P, C, S> = CompressedShell<P, C, S>;
+
+trait SP<S>: Fn(&S, Point3, Option<(f64, f64)>) -> Option<(f64, f64)> {}
+impl<S, F> SP<S> for F where F: Fn(&S, Point3, Option<(f64, f64)>) -> Option<(f64, f64)> {}
+
+mod split_closed_edges;
+use split_closed_edges::split_closed_edges;
+
+mod split_closed_faces;
+use split_closed_faces::split_closed_faces;
+
+/// Splits closed edges and faces
+///
+/// # Details
+/// The topology of the shapes handled by truck has the following rules
+/// - The endpoints of the edges must be different.
+/// - The boundaries of the faces must be a simple wire.
+///
+/// Shapes created in other CAD systems do not necessarily follow these rules.
+/// When such shapes are handled by truck, this method is applied at the stage
+/// of `CompressedShell` and `CompressedSolid`, which are intermediate forms.
+///
+/// # Remarks
+/// Boundary simplification is still only implemented for cylinders.
+/// It has not yet been implemented for cases involving singularities, such as spherical surfaces.
+pub trait SplitClosedEdgesAndFaces {
+    /// Splits closed edges and faces
+    fn split_closed_edges_and_faces(&mut self, tol: f64);
+}
+
+impl<C, S> SplitClosedEdgesAndFaces for CompressedShell<Point3, C, S>
+where
+    C: ParametricCurve3D
+        + BoundedCurve
+        + Cut
+        + ParameterDivision1D<Point = Point3>
+        + SearchNearestParameter<D1, Point = Point3>
+        + TryFrom<PCurve<Line<Point2>, S>>,
+    S: ParametricSurface3D + SearchParameter<D2, Point = Point3>,
+{
+    fn split_closed_edges_and_faces(&mut self, tol: f64) {
+        fn sp<S>(surface: &S, point: Point3, hint: Option<(f64, f64)>) -> Option<(f64, f64)>
+        where S: SearchParameter<D2, Point = Point3> {
+            surface.search_parameter(point, hint, 100)
+        }
+        split_closed_edges(self);
+        split_closed_faces(self, tol, sp);
+    }
+}
+
+impl<C, S> SplitClosedEdgesAndFaces for CompressedSolid<Point3, C, S>
+where
+    C: ParametricCurve3D
+        + BoundedCurve
+        + Cut
+        + ParameterDivision1D<Point = Point3>
+        + SearchNearestParameter<D1, Point = Point3>
+        + TryFrom<PCurve<Line<Point2>, S>>,
+    S: ParametricSurface3D + SearchParameter<D2, Point = Point3>,
+{
+    fn split_closed_edges_and_faces(&mut self, tol: f64) {
+        self.boundaries
+            .iter_mut()
+            .for_each(|shell| shell.split_closed_edges_and_faces(tol))
+    }
+}
+
+/// robust version of splitting closed edges and faces.
+///
+/// # Details
+/// Robust version of [`SplitClosedEdgesAndFaces`] based on [`SearchNearestParameter`].
+pub trait RobustSplitClosedEdgesAndFaces {
+    /// Splits closed edges and faces
+    fn robust_split_closed_edges_and_faces(&mut self, tol: f64);
+}
+
+impl<C, S> RobustSplitClosedEdgesAndFaces for CompressedShell<Point3, C, S>
+where
+    C: ParametricCurve3D
+        + BoundedCurve
+        + Cut
+        + ParameterDivision1D<Point = Point3>
+        + SearchNearestParameter<D1, Point = Point3>
+        + TryFrom<PCurve<Line<Point2>, S>>,
+    S: ParametricSurface3D
+        + SearchParameter<D2, Point = Point3>
+        + SearchNearestParameter<D2, Point = Point3>,
+{
+    fn robust_split_closed_edges_and_faces(&mut self, tol: f64) {
+        fn sp<S>(surface: &S, point: Point3, hint: Option<(f64, f64)>) -> Option<(f64, f64)>
+        where S: SearchParameter<D2, Point = Point3> + SearchNearestParameter<D2, Point = Point3>
+        {
+            surface
+                .search_parameter(point, hint, 100)
+                .or_else(|| surface.search_parameter(point, None, 100))
+                .or_else(|| surface.search_nearest_parameter(point, hint, 100))
+                .or_else(|| surface.search_nearest_parameter(point, None, 100))
+        }
+        split_closed_edges(self);
+        split_closed_faces(self, tol, sp);
+    }
+}
+
+impl<C, S> RobustSplitClosedEdgesAndFaces for CompressedSolid<Point3, C, S>
+where
+    C: ParametricCurve3D
+        + BoundedCurve
+        + Cut
+        + ParameterDivision1D<Point = Point3>
+        + SearchNearestParameter<D1, Point = Point3>
+        + TryFrom<PCurve<Line<Point2>, S>>,
+    S: ParametricSurface3D
+        + SearchParameter<D2, Point = Point3>
+        + SearchNearestParameter<D2, Point = Point3>,
+{
+    fn robust_split_closed_edges_and_faces(&mut self, tol: f64) {
+        let fs = RobustSplitClosedEdgesAndFaces::robust_split_closed_edges_and_faces;
+        self.boundaries.iter_mut().for_each(|shell| fs(shell, tol))
+    }
+}
+
+#[cfg(test)]
+mod tests;
@@ -0,0 +1,68 @@
+use super::*;
+
+/// Cuts the edge at the middle point of the bounded curve `curve`.
+fn cut_edge<P, C>(
+    vertices: &mut Vec<P>,
+    Edge {
+        vertices: (_, v1),
+        curve,
+    }: &mut Edge<C>,
+) -> Edge<C>
+where
+    C: BoundedCurve<Point = P> + Cut,
+{
+    let (t0, t1) = curve.range_tuple();
+    let t = (t0 + t1) / 2.0;
+    vertices.push(curve.subs(t));
+
+    let v2 = *v1;
+    *v1 = vertices.len() - 1;
+    let v1 = *v1;
+    let curve0 = curve.cut(t);
+    Edge {
+        vertices: (v1, v2),
+        curve: curve0,
+    }
+}
+
+/// Adds cut edges to the vector `edges`.
+fn add_edges<P, C>(vertices: &mut Vec<P>, edges: &mut Vec<Edge<C>>) -> HashMap<usize, usize>
+where C: BoundedCurve<Point = P> + Cut {
+    let len = edges.len();
+    let sub_add_edges = move |i: usize| {
+        let (v0, v1) = edges[i].vertices;
+        if v0 == v1 {
+            let new_edge = cut_edge(vertices, &mut edges[i]);
+            edges.push(new_edge);
+            Some((i, edges.len() - 1))
+        } else {
+            None
+        }
+    };
+    (0..len).filter_map(sub_add_edges).collect()
+}
+
+fn replace_edges(wire: &mut Wire, added: &HashMap<usize, usize>) {
+    let insert_one_edge = |edge: EdgeIndex| {
+        let Some(new_edge_index) = added.get(&edge.index) else {
+            return vec![edge];
+        };
+        let new_edge = EdgeIndex {
+            index: *new_edge_index,
+            orientation: edge.orientation,
+        };
+        match edge.orientation {
+            true => vec![edge, new_edge],
+            false => vec![new_edge, edge],
+        }
+    };
+    let new_wire = wire.iter().copied().flat_map(insert_one_edge).collect();
+    *wire = new_wire;
+}
+
+pub(super) fn split_closed_edges<P, C, S>(shell: &mut Shell<P, C, S>)
+where C: BoundedCurve<Point = P> + Cut {
+    let added = add_edges(&mut shell.vertices, &mut shell.edges);
+    let wire_iter = shell.faces.iter_mut().flat_map(|face| &mut face.boundaries);
+    wire_iter.for_each(|wire| replace_edges(wire, &added));
+}
@@ -0,0 +1,975 @@
+use super::*;
+use std::ops::Range;
+
+pub(super) fn split_closed_faces<C, S>(shell: &mut Shell<Point3, C, S>, tol: f64, sp: impl SP<S>)
+where
+    C: ParametricCurve3D
+        + BoundedCurve
+        + Cut
+        + ParameterDivision1D<Point = Point3>
+        + SearchNearestParameter<D1, Point = Point3>
+        + TryFrom<PCurve<Line<Point2>, S>>,
+    S: ParametricSurface3D, {
+    let to_poly = closure_to_poly(tol);
+    let mut poly_edges: Vec<_> = shell.edges.iter().map(to_poly).collect();
+    let len = shell.faces.len();
+    (0..len).for_each(|i| {
+        split_face_with_non_closed_boundary(i, shell, &mut poly_edges, &sp, tol);
+    });
+    let len = shell.faces.len();
+    let closure = |i| split_face_with_non_simple_wire(i, shell, &mut poly_edges, &sp, tol);
+    let new_faces: Vec<_> = (0..len).filter_map(closure).flatten().collect();
+    shell.faces.extend(new_faces);
+}
+
+fn split_face_with_non_closed_boundary<C, S>(
+    face_index: usize,
+    shell: &mut Shell<Point3, C, S>,
+    poly_edges: &mut Vec<PolylineCurve<Point3>>,
+    sp: impl SP<S>,
+    tol: f64,
+) -> Option<()>
+where
+    C: ParametricCurve3D
+        + BoundedCurve
+        + Cut
+        + ParameterDivision1D<Point = Point3>
+        + SearchNearestParameter<D1, Point = Point3>
+        + TryFrom<PCurve<Line<Point2>, S>>,
+    S: ParametricSurface3D,
+{
+    let (divisor, open, closed) =
+        non_closed_wires_in_param_divisor(face_index, shell, poly_edges, &sp)?;
+    debug_assert_eq!(open.len(), 2);
+    take_vertices_to_intersections(divisor, face_index, shell, poly_edges, &sp, tol);
+    let Shell { faces, edges, .. } = shell;
+    let new_boundaries = split_boundaries_by_divisor(
+        &faces[face_index],
+        &closed,
+        divisor,
+        edges,
+        poly_edges,
+        &sp,
+        tol,
+    )?;
+    let boundaries = &mut faces[face_index].boundaries;
+    connect_open_boundaries(boundaries, new_boundaries, divisor, edges, open)
+}
+
+fn split_face_with_non_simple_wire<C, S>(
+    face_index: usize,
+    shell: &mut Shell<Point3, C, S>,
+    poly_edges: &mut Vec<PolylineCurve<Point3>>,
+    sp: impl SP<S>,
+    tol: f64,
+) -> Option<Vec<Face<S>>>
+where
+    C: ParametricCurve3D
+        + BoundedCurve
+        + Cut
+        + ParameterDivision1D<Point = Point3>
+        + SearchNearestParameter<D1, Point = Point3>
+        + TryFrom<PCurve<Line<Point2>, S>>,
+    S: ParametricSurface3D,
+{
+    let Face { boundaries, .. } = &shell.faces[face_index];
+    let divisor = find_non_simple_wire_divisor(boundaries, &shell.edges, &shell.vertices)?;
+    split_face_by_divisor(face_index, divisor, shell, poly_edges, sp, tol)
+}
+
+// --- find_non_simple_wire_divisor ---
+
+fn find_non_simple_wire_divisor<C>(
+    boundaries: &[Wire],
+    edges: &[Edge<C>],
+    vertices: &[Point3],
+) -> Option<(usize, usize)> {
+    let closure = |boundary| non_simple_wire_divisor(boundary, edges, vertices);
+    boundaries.iter().cloned().find_map(closure)
+}
+
+fn non_simple_wire_divisor<C>(
+    mut boundary: Wire,
+    edges: &[Edge<C>],
+    vertices: &[Point3],
+) -> Option<(usize, usize)> {
+    let (i0, j0) = boundary.iter().enumerate().find_map(find_loop())?;
+    boundary.rotate_left(i0);
+    let j0 = j0 - i0;
+
+    let (i1, j1) = boundary
+        .iter()
+        .enumerate()
+        .skip(j0 + 1)
+        .find_map(find_loop())
+        .unwrap_or((j0, boundary.len()));
+
+    let (k0, k1) = ((j0 + 1) / 2, (i1 + j1 + 1) / 2);
+    let f = closure_take_front(edges);
+    let pre_divisor = (f(boundary[k0]), f(boundary[k1]));
+    nearest_correction(
+        pre_divisor,
+        &boundary[1..j0],
+        &boundary[i1 + 1..j1],
+        vertices,
+        closure_take_front(edges),
+    )
+}
+
+fn find_loop() -> impl FnMut((usize, &EdgeIndex)) -> Option<(usize, usize)> {
+    let mut map = HashMap::<usize, usize>::default();
+    move |(i, edge)| map.insert(edge.index, i).map(move |j| (j, i))
+}
+
+// --- find_non_closed_wires_in_param_divisor ---
+
+type FindNonClosedWiresInParamDivisorResult = Option<((usize, usize), Vec<Wire>, Vec<Wire>)>;
+fn non_closed_wires_in_param_divisor<C, S>(
+    face_index: usize,
+    Shell {
+        edges,
+        faces,
+        ref vertices,
+    }: &mut Shell<Point3, C, S>,
+    poly_edges: &[PolylineCurve<Point3>],
+    sp: impl SP<S>,
+) -> FindNonClosedWiresInParamDivisorResult
+where
+    C: ParametricCurve3D + BoundedCurve + TryFrom<PCurve<Line<Point2>, S>>,
+    S: ParametricSurface3D,
+{
+    let Face {
+        boundaries,
+        surface,
+        ..
+    } = &faces[face_index];
+    let param_boundaries = create_param_boundaries(boundaries, surface, poly_edges, sp)?;
+    let (mut open, mut closed) = (Vec::new(), Vec::new());
+    boundaries
+        .iter()
+        .zip(&param_boundaries)
+        .for_each(|(boundary, param_boundary)| {
+            let poly = PolylineCurve(param_boundary.iter().flatten().copied().collect());
+            match boundary_type(&poly) {
+                BoundaryType::NotClosed => open.push((boundary, param_boundary)),
+                _ => closed.push((boundary, param_boundary)),
+            }
+        });
+    if open.len() != 2 {
+        return None;
+    }
+    let take_front = &closure_take_front(edges);
+    let find_nearest_to_ends = move |i: usize, u0: f64, up: f64| {
+        move |(edge_index, poly): ZippedEdge<'_>| {
+            let u = poly[0][i];
+            let closure = move |i| {
+                let u = u + i as f64 * up;
+                f64::abs(u - u0)
+            };
+            let dist = (-2..=2)
+                .map(closure)
+                .min_by(|x, y| x.partial_cmp(y).unwrap())
+                .unwrap();
+            (take_front(*edge_index), dist)
+        }
+    };
+    let create_pre_divisor = move |(wire, poly_wire): ZippedWire<'_>| {
+        let (urange, vrange) = surface.try_range_tuple();
+        let (up, vp) = (surface.u_period(), surface.v_period());
+        let p = poly_wire[0][0];
+        let q = *poly_wire.last().unwrap().last().unwrap();
+        let vertices: Vec<(usize, f64)> = if p.x.near(&q.x) {
+            if let (Some(vp), Some((v0, _))) = (vp, vrange) {
+                let closure = find_nearest_to_ends(1, v0, vp);
+                Some(wire.iter().zip(poly_wire).map(closure).collect())
+            } else {
+                None
+            }
+        } else if p.y.near(&q.y) {
+            if let (Some(up), Some((u0, _))) = (up, urange) {
+                let closure = find_nearest_to_ends(0, u0, up);
+                Some(wire.iter().zip(poly_wire).map(closure).collect())
+            } else {
+                None
+            }
+        } else {
+            None
+        }?;
+        vertices
+            .into_iter()
+            .min_by(|(_, x), (_, y)| x.partial_cmp(y).unwrap())
+            .map(|(v, _)| v)
+    };
+    let pre_divisor_vec = open
+        .iter()
+        .copied()
+        .map(create_pre_divisor)
+        .collect::<Option<Vec<_>>>()?;
+    let open: Vec<_> = open.into_iter().map(|(w, _)| w.to_vec()).collect();
+    let closed: Vec<_> = closed.into_iter().map(|(w, _)| w.to_vec()).collect();
+    let divisor = nearest_correction(
+        (pre_divisor_vec[0], pre_divisor_vec[1]),
+        &open[0],
+        &open[1],
+        vertices,
+        take_front,
+    )?;
+    Some((divisor, open, closed))
+}
+
+// --- split_face_by_divisor ---
+
+fn split_face_by_divisor<C, S>(
+    face_index: usize,
+    divisor: (usize, usize),
+    shell: &mut Shell<Point3, C, S>,
+    poly_edges: &mut Vec<PolylineCurve<Point3>>,
+    sp: impl SP<S>,
+    tol: f64,
+) -> Option<Vec<Face<S>>>
+where
+    C: ParametricCurve3D
+        + BoundedCurve
+        + Cut
+        + ParameterDivision1D<Point = Point3>
+        + SearchNearestParameter<D1, Point = Point3>
+        + TryFrom<PCurve<Line<Point2>, S>>,
+    S: ParametricSurface3D,
+{
+    take_vertices_to_intersections(divisor, face_index, shell, poly_edges, &sp, tol);
+    let face = &shell.faces[face_index];
+    let new_boundaries = split_boundaries_by_divisor(
+        face,
+        &face.boundaries,
+        divisor,
+        &mut shell.edges,
+        poly_edges,
+        &sp,
+        tol,
+    )?;
+    divide_face(
+        &mut shell.faces[face_index],
+        new_boundaries,
+        poly_edges,
+        &sp,
+    )
+}
+
+// --- take_vertices_to_intersections ---
+
+fn take_vertices_to_intersections<C, S>(
+    divisor: (usize, usize),
+    face_index: usize,
+    Shell {
+        vertices,
+        edges,
+        faces,
+    }: &mut Shell<Point3, C, S>,
+    poly_edges: &mut Vec<PolylineCurve<Point3>>,
+    sp: impl SP<S>,
+    tol: f64,
+) -> Option<()>
+where
+    C: ParametricCurve3D
+        + BoundedCurve
+        + Cut
+        + SearchNearestParameter<D1, Point = Point3>
+        + ParameterDivision1D<Point = Point3>,
+    S: ParametricSurface3D,
+{
+    let Face {
+        ref boundaries,
+        ref surface,
+        ..
+    } = faces[face_index];
+    let param_boundaries = create_param_boundaries(boundaries, surface, poly_edges, &sp)?;
+    let (v0, v1) = divisor;
+    let p = zip_boundaries(boundaries, &param_boundaries)
+        .find_map(closure_find_vertex_parameter(v0, edges))?;
+    let q = zip_boundaries(boundaries, &param_boundaries)
+        .find_map(closure_find_vertex_parameter(v1, edges))?;
+    let Face {
+        ref boundaries,
+        ref surface,
+        ..
+    } = faces[face_index];
+    let periods = (surface.u_period(), surface.v_period());
+    let q = periodic_iterator(q, periods)
+        .min_by(|q, r| p.distance2(*q).partial_cmp(&p.distance2(*r)).unwrap())?;
+    let pcurve = PCurve::new(Line(p, q), surface.clone());
+    let cut_edge = |(edge_index, param_edge): ZippedEdge<'_>| {
+        let index = edge_index.index;
+        let vec = enumerate_intersections(&edges[index], param_edge, &pcurve)?;
+        if vec.is_empty() {
+            return Some(None);
+        }
+        let erange = cut_edge_by_intersections(index, vec, edges, poly_edges, vertices, tol);
+        Some(Some((index, erange)))
+    };
+    let new_edges = zip_boundaries(boundaries, &param_boundaries)
+        .map(cut_edge)
+        .collect::<Option<Vec<_>>>()?;
+    let insert = |(index, erange): (usize, Range<usize>)| {
+        faces
+            .iter_mut()
+            .flat_map(|face| &mut face.boundaries)
+            .for_each(|wire| insert_new_edges(wire, index, erange.clone()));
+    };
+    new_edges.into_iter().flatten().for_each(insert);
+    Some(())
+}
+
+fn closure_find_vertex_parameter<C>(
+    v: usize,
+    edges: &[Edge<C>],
+) -> impl Fn(ZippedEdge<'_>) -> Option<Point2> + '_ {
+    let take_front = closure_take_front(edges);
+    move |(edge_index, param_edge): ZippedEdge<'_>| match v == take_front(*edge_index) {
+        true => Some(param_edge[0]),
+        false => None,
+    }
+}
+
+fn enumerate_intersections<C, S>(
+    edge: &Edge<C>,
+    param_edge: &PolylineCurve<Point2>,
+    pcurve: &PCurve<Line<Point2>, S>,
+) -> Option<Vec<(f64, Point3)>>
+where
+    C: ParametricCurve3D + SearchNearestParameter<D1, Point = Point3>,
+    S: ParametricSurface3D,
+{
+    let rough_intersections = intersections_between_line_polyline(*pcurve.curve(), param_edge);
+    let intersections = exact_intersections(rough_intersections, pcurve, &edge.curve)?;
+    Some(intersections)
+}
+
+fn intersections_between_line_polyline(
+    line: Line<Point2>,
+    param_edge: &PolylineCurve<Point2>,
+) -> Vec<f64> {
+    let filter = |p: &[Point2]| {
+        let (s, t, _) = line.intersection(Line(p[0], p[1]))?;
+        let unit = 0.0..1.0;
+        match unit.contains(&s) && unit.contains(&t) {
+            true => Some(s),
+            false => None,
+        }
+    };
+    let mut vec: Vec<_> = param_edge.windows(2).filter_map(filter).collect();
+    if !vec.is_empty() {
+        let first = *param_edge.first().unwrap();
+        let last = *param_edge.last().unwrap();
+        if line.distance_to_point(first).so_small() {
+            vec.remove(0);
+        }
+        if line.distance_to_point(last).so_small() {
+            vec.pop();
+        }
+    }
+    vec
+}
+
+fn exact_intersections<C, S>(
+    naive_intersections: Vec<f64>,
+    pcurve: &PCurve<Line<Point2>, S>,
+    curve: &C,
+) -> Option<Vec<(f64, Point3)>>
+where
+    C: ParametricCurve3D + SearchNearestParameter<D1, Point = Point3>,
+    S: ParametricSurface3D,
+{
+    let closure = |t| {
+        let (_, t, p) = search_intersection(pcurve, curve, t)?;
+        Some((t, p))
+    };
+    let mut res = naive_intersections
+        .into_iter()
+        .map(closure)
+        .collect::<Option<Vec<_>>>()?;
+    res.sort_by(|(s, _), (t, _)| s.partial_cmp(t).unwrap());
+    Some(res)
+}
+
+fn search_intersection<C, S>(
+    pcurve: &PCurve<Line<Point2>, S>,
+    curve: &C,
+    mut t0: f64,
+) -> Option<(f64, f64, Point3)>
+where
+    C: ParametricCurve3D + SearchNearestParameter<D1, Point = Point3>,
+    S: ParametricSurface3D,
+{
+    let mut previous0 = t0;
+    let mut previous1 = None;
+    for _ in 0..100 {
+        let p = pcurve.subs(t0);
+        let t1 = curve.search_nearest_parameter(p, previous1, 100)?;
+        let q = curve.subs(t1);
+        t0 = pcurve.search_nearest_parameter(q, t0, 100)?;
+        if let Some(previous1) = previous1 {
+            if previous0.near(&t0) && previous1.near(&t1) {
+                return Some((t0, t1, pcurve.subs(t0)));
+            }
+        }
+        previous0 = t0;
+        previous1 = Some(t1);
+    }
+    None
+}
+
+fn cut_edge_by_intersections<C>(
+    edge_index: usize,
+    intersections: Vec<(f64, Point3)>,
+    edges: &mut Vec<Edge<C>>,
+    poly_edges: &mut Vec<PolylineCurve<Point3>>,
+    vertices: &mut Vec<Point3>,
+    tol: f64,
+) -> Range<usize>
+where
+    C: Cut<Point = Point3> + ParameterDivision1D<Point = Point3>,
+{
+    let edge = &mut edges[edge_index];
+    let vfirst = vertices.len();
+    vertices.extend(intersections.iter().map(|(_, p)| *p));
+    let intersections_iter = intersections.into_iter().enumerate().rev();
+    let closure = |(i, (t, _))| Edge {
+        vertices: (vfirst + i, vfirst + i + 1),
+        curve: edge.curve.cut(t),
+    };
+    let mut new_edges: Vec<_> = intersections_iter.map(closure).collect();
+    let to_poly = closure_to_poly(tol);
+    poly_edges[edge_index] = to_poly(edge);
+    poly_edges.extend(new_edges.iter().map(to_poly));
+    new_edges.last_mut().unwrap().vertices.1 = edge.vertices.1;
+    edge.vertices.1 = vfirst;
+    let efirst = edges.len();
+    edges.extend(new_edges);
+    efirst..edges.len()
+}
+
+fn insert_new_edges(wire: &mut Wire, pivot_edge_index: usize, inserted_range: Range<usize>) {
+    let positions = wire
+        .iter()
+        .enumerate()
+        .filter_map(
+            |(i, edge_index)| match edge_index.index == pivot_edge_index {
+                true => Some(i),
+                false => None,
+            },
+        )
+        .collect::<Vec<usize>>();
+    positions.into_iter().rev().for_each(|i| {
+        let orientation = wire[i].orientation;
+        let ei = |index| EdgeIndex { index, orientation };
+        if orientation {
+            let tmp_wire = wire.split_off(i + 1);
+            wire.extend(inserted_range.clone().map(ei));
+            wire.extend(tmp_wire);
+        } else {
+            let tmp_wire = wire.split_off(i);
+            wire.extend(inserted_range.clone().rev().map(ei));
+            wire.extend(tmp_wire);
+        }
+    });
+}
+
+// --- split boundaries ---
+
+fn split_boundaries_by_divisor<C, S>(
+    Face {
+        boundaries,
+        surface,
+        ..
+    }: &Face<S>,
+    closed: &[Wire],
+    divisor: (usize, usize),
+    edges: &mut Vec<Edge<C>>,
+    poly_edges: &mut Vec<PolylineCurve<Point3>>,
+    sp: impl SP<S>,
+    tol: f64,
+) -> Option<Vec<Wire>>
+where
+    C: ParametricCurve3D + BoundedCurve + TryFrom<PCurve<Line<Point2>, S>>,
+    S: ParametricSurface3D,
+{
+    let param_boundaries = create_param_boundaries(boundaries, surface, poly_edges, &sp)?;
+    let param_vertices = create_param_vertices(boundaries, &param_boundaries, edges);
+    let mut duplicated_edges = duplicated_edges(boundaries.iter().flatten().map(|ei| ei.index));
+    let vertices_on_divisor = enumerate_vertices_on_divisor(divisor, &param_vertices, surface)?;
+    let new_edges = create_new_edges(&vertices_on_divisor, poly_edges, surface, tol)?;
+    let (new_edge_range, new_edge_indices) = signup_new_edges(edges, new_edges);
+    duplicated_edges.extend(new_edge_range);
+    let edge_iter = closed.iter().flatten().copied().chain(new_edge_indices);
+    let vemap = create_vemap(edges, edge_iter);
+    Some(construct_boundaries(vemap, edges, &duplicated_edges))
+}
+
+// --- create_param_boundaries ---
+
+fn create_param_boundaries<S: ParametricSurface3D>(
+    boundaries: &[Wire],
+    surface: &S,
+    poly_edges: &[PolylineCurve<Point3>],
+    sp: impl SP<S>,
+) -> Option<Vec<Vec<PolylineCurve<Point2>>>> {
+    let closure = |boundary: &Wire| boundary_to_param_polys(boundary, surface, poly_edges, &sp);
+    boundaries.iter().map(closure).collect()
+}
+
+fn boundary_to_param_polys<S: ParametricSurface3D>(
+    boundary: &Wire,
+    surface: &S,
+    poly_edges: &[PolylineCurve<Point3>],
+    sp: impl SP<S>,
+) -> Option<Vec<PolylineCurve<Point2>>> {
+    let get_poly = closure_get_poly(poly_edges);
+    let poly_boundary: Vec<_> = boundary.iter().copied().map(get_poly).collect();
+    let poly_wire_iter = PolyWireIter::try_new(&poly_boundary)?;
+    let closure = poly_project_to_uv(surface, sp);
+    let mut long_poly = poly_wire_iter.map(closure).collect::<Option<Vec<_>>>()?;
+    boundary_into_domain(&mut long_poly, surface);
+    let mut vec = vec![long_poly.remove(0)];
+    let split_into_2dpoly = move |poly: &PolylineCurve<_>| {
+        let latter = long_poly.split_off(poly.len() - 1);
+        vec.append(&mut long_poly);
+        let mut res = Vec::new();
+        res.append(&mut vec);
+        vec.push(*res.last().unwrap());
+        long_poly = latter;
+        PolylineCurve(res)
+    };
+    Some(poly_boundary.iter().map(split_into_2dpoly).collect())
+}
+
+struct PolyWireIter<'a, P> {
+    current: std::slice::Iter<'a, P>,
+    stock: std::vec::IntoIter<std::slice::Iter<'a, P>>,
+}
+
+impl<'a, P> PolyWireIter<'a, P> {
+    fn try_new<I: AsRef<[P]> + 'a>(vec: &'a [I]) -> Option<Self> {
+        let iters: Vec<_> = vec.iter().map(|v| v.as_ref().iter()).collect();
+        let mut stock = iters.into_iter();
+        Some(Self {
+            current: stock.next()?,
+            stock,
+        })
+    }
+}
+
+impl<'a, P: Copy> Iterator for PolyWireIter<'a, P> {
+    type Item = P;
+    fn next(&mut self) -> Option<Self::Item> {
+        if let Some(next) = self.current.next() {
+            Some(*next)
+        } else {
+            self.current = self.stock.next()?;
+            self.current.next()?;
+            self.current.next().cloned()
+        }
+    }
+}
+
+fn poly_project_to_uv<'a, S: ParametricSurface3D>(
+    surface: &'a S,
+    sp: impl SP<S> + 'a,
+) -> impl FnMut(Point3) -> Option<Point2> + 'a {
+    let (up, vp) = (surface.u_period(), surface.v_period());
+    let mut previous = None;
+    move |pt| {
+        let (mut u, mut v) = sp(surface, pt, previous)?;
+        if let (Some(up), Some((u0, _))) = (up, previous) {
+            u = get_mindiff(u, u0, up);
+        }
+        if let (Some(vp), Some((_, v0))) = (vp, previous) {
+            v = get_mindiff(v, v0, vp);
+        }
+        previous = Some((u, v));
+        Some(Point2::new(u, v))
+    }
+}
+
+fn abs_diff(previous: f64) -> impl Fn(&f64, &f64) -> std::cmp::Ordering {
+    let f = move |x: &f64| f64::abs(x - previous);
+    move |x: &f64, y: &f64| f(x).partial_cmp(&f(y)).unwrap()
+}
+fn get_mindiff(u: f64, u0: f64, up: f64) -> f64 {
+    let closure = |i| u + i as f64 * up;
+    (-2..=2).map(closure).min_by(abs_diff(u0)).unwrap()
+}
+
+fn boundary_into_domain<S: ParametricSurface3D>(vec: &mut Vec<Point2>, surface: &S) {
+    let (up, vp) = (surface.u_period(), surface.v_period());
+    let (urange, vrange) = surface.try_range_tuple();
+    let grav = vec.iter().fold(Point2::origin(), |g, p| g + p.to_vec()) / vec.len() as f64;
+    if let (Some(up), Some((u0, _))) = (up, urange) {
+        let quot = f64::floor((grav.x - u0) / up);
+        vec.iter_mut().for_each(|p| p.x -= quot * up);
+    }
+    if let (Some(vp), Some((v0, _))) = (vp, vrange) {
+        let quot = f64::floor((grav.y - v0) / vp);
+        vec.iter_mut().for_each(|p| p.y -= quot * vp);
+    }
+    let last = *vec.last().unwrap();
+    if !vec[0].near(&last) {
+        let Point2 { x: u0, y: v0 } = last;
+        if surface.uder(u0, v0).so_small() || surface.vder(u0, v0).so_small() {
+            vec.push(vec[0]);
+        }
+    }
+}
+
+// --- create_param_vertices ---
+
+fn create_param_vertices<C>(
+    boundaries: &[Wire],
+    param_boundaries: &[Vec<PolylineCurve<Point2>>],
+    edges: &[Edge<C>],
+) -> HashMap<usize, Point2> {
+    let take_front = closure_take_front(edges);
+    zip_boundaries(boundaries, param_boundaries)
+        .map(|(edge_index, param_edge)| (take_front(*edge_index), param_edge[0]))
+        .collect()
+}
+
+// --- duplicated_edges ---
+
+fn duplicated_edges(edges: impl Iterator<Item = usize>) -> HashSet<usize> {
+    let mut done = HashSet::default();
+    edges.filter(move |index| !done.insert(*index)).collect()
+}
+
+// --- enumerate_vertices_on_divisor ---
+
+fn enumerate_vertices_on_divisor<S: ParametricSurface>(
+    divisor: (usize, usize),
+    param_vertices: &HashMap<usize, Point2>,
+    surface: &S,
+) -> Option<Vec<(f64, (usize, Point2))>> {
+    let (v0, v1) = divisor;
+    let (p, q) = (*param_vertices.get(&v0)?, *param_vertices.get(&v1)?);
+    let periods = (surface.u_period(), surface.v_period());
+    let q = periodic_iterator(q, periods)
+        .min_by(|q, r| p.distance2(*q).partial_cmp(&p.distance2(*r)).unwrap())?;
+    let line = Line(p, q);
+    let iter = param_vertices.iter().filter_map(move |(v, uv)| {
+        periodic_iterator(*uv, periods)
+            .find(move |uv| line.distance_to_point_as_segment(*uv).so_small())
+            .map(move |uv| Some((line.search_nearest_parameter(uv, None, 1)?, (*v, uv))))
+    });
+    let mut vertices_on_divisor = iter.collect::<Option<Vec<_>>>()?;
+    vertices_on_divisor.sort_by(|(s, _), (t, _)| s.partial_cmp(t).unwrap());
+    Some(vertices_on_divisor)
+}
+
+// --- create_new_edges ---
+
+fn create_new_edges<C, S>(
+    vertices_on_divisor: &[(f64, (usize, Point2))],
+    poly_edges: &mut Vec<PolylineCurve<Point3>>,
+    surface: &S,
+    tol: f64,
+) -> Option<Vec<Edge<C>>>
+where
+    C: TryFrom<PCurve<Line<Point2>, S>>,
+    S: ParametricSurface3D,
+{
+    vertices_on_divisor
+        .iter()
+        .for_each(|(_, (_, uv))| print!("{:?} ", surface.subs(uv.x, uv.y)));
+    println!();
+    let make_edge = move |p: &[(f64, (usize, Point2))]| {
+        let ((_, (v0, uv0)), (_, (v1, uv1))) = (p[0], p[1]);
+        let pcurve = PCurve::new(Line(uv0, uv1), surface.clone());
+        poly_edges.push(PolylineCurve::from_curve(&pcurve, (0.0, 1.0), tol));
+        Some(Edge {
+            vertices: (v0, v1),
+            curve: C::try_from(pcurve).ok()?,
+        })
+    };
+    vertices_on_divisor.chunks(2).map(make_edge).collect()
+}
+
+// --- signup new edges ---
+
+fn signup_new_edges<C>(
+    edges: &mut Vec<Edge<C>>,
+    new_edges: Vec<Edge<C>>,
+) -> (Range<usize>, impl Iterator<Item = EdgeIndex>) {
+    let len = edges.len();
+    edges.extend(new_edges);
+    let new_edge_indices = (len..edges.len()).flat_map(move |index| {
+        let ei = |orientation: bool| EdgeIndex { index, orientation };
+        [ei(true), ei(false)]
+    });
+    (len..edges.len(), new_edge_indices)
+}
+
+// --- create_vemap ---
+
+fn create_vemap<C>(
+    edges: &[Edge<C>],
+    edge_iter: impl Iterator<Item = EdgeIndex>,
+) -> HashMap<usize, Vec<EdgeIndex>> {
+    let mut vemap = HashMap::<usize, Vec<EdgeIndex>>::default();
+    let take_front = closure_take_front(edges);
+    edge_iter.for_each(|edge_index| {
+        let v = take_front(edge_index);
+        vemap.entry(v).or_insert_with(Vec::new).push(edge_index);
+    });
+    vemap
+}
+
+// --- construct boundaries ---
+
+fn construct_boundaries<C>(
+    mut vemap: HashMap<usize, Vec<EdgeIndex>>,
+    edges: &[Edge<C>],
+    new_edge_range: &HashSet<usize>,
+) -> Vec<Wire> {
+    let take_back = closure_take_back(edges);
+    let mut new_boundaries = Vec::new();
+    while !vemap.is_empty() {
+        let (start, vec) = vemap.iter_mut().min_by_key(|(idx, _)| *idx).unwrap();
+        let start = *start;
+        let mut edge_index = vec.pop().unwrap();
+        if vec.is_empty() {
+            vemap.remove(&start);
+        }
+        let mut wire = Vec::new();
+        loop {
+            wire.push(edge_index);
+            let v = take_back(edge_index);
+            if v == start {
+                break;
+            }
+            let Some(vec) = vemap.get_mut(&v) else {
+                unreachable!();
+            };
+            match vec.len() {
+                1 => {
+                    edge_index = vec.pop().unwrap();
+                    vemap.remove(&v);
+                }
+                2 => {
+                    let is_new0 = new_edge_range.contains(&edge_index.index);
+                    let is_new1 = new_edge_range.contains(&vec[0].index);
+                    let i = if is_new0 != is_new1 { 0 } else { 1 };
+                    edge_index = vec.remove(i);
+                }
+                _ => panic!("something wrong!"),
+            }
+        }
+        new_boundaries.push(wire);
+    }
+    new_boundaries
+}
+
+// --- divide_face ---
+
+fn divide_face<S: ParametricSurface3D>(
+    Face {
+        ref mut boundaries,
+        ref surface,
+        ref orientation,
+    }: &mut Face<S>,
+    new_boundaries: Vec<Wire>,
+    poly_edges: &[PolylineCurve<Point3>],
+    sp: impl SP<S>,
+) -> Option<Vec<Face<S>>> {
+    let mut face_boundaries = assort_boundary(surface, new_boundaries, poly_edges, sp)?.into_iter();
+    *boundaries = face_boundaries.next().unwrap();
+    let create_face = |boundaries| Face {
+        boundaries,
+        surface: surface.clone(),
+        orientation: *orientation,
+    };
+    Some(face_boundaries.map(create_face).collect())
+}
+
+fn assort_boundary<S>(
+    surface: &S,
+    boundaries: Vec<Wire>,
+    poly_edges: &[PolylineCurve<Point3>],
+    sp: impl SP<S>,
+) -> Option<Vec<Vec<Wire>>>
+where
+    S: ParametricSurface3D,
+{
+    let (mut positives, mut negatives) = (Vec::new(), Vec::new());
+    let param_boundaries = create_param_boundaries(&boundaries, surface, poly_edges, sp)?;
+    boundaries
+        .into_iter()
+        .zip(param_boundaries)
+        .try_for_each(|(boundary, param_boundary)| {
+            let poly_boundary = PolylineCurve(param_boundary.into_iter().flatten().collect());
+            match boundary_type(&poly_boundary) {
+                BoundaryType::Positive => positives.push((boundary, poly_boundary)),
+                BoundaryType::Negative => negatives.push((boundary, poly_boundary)),
+                _ => return None,
+            }
+            Some(())
+        })?;
+    let mut assorted = positives.into_iter().map(|x| vec![x]).collect::<Vec<_>>();
+    negatives
+        .into_iter()
+        .try_for_each(|(boundary, poly_boundary)| {
+            let p = poly_boundary[0];
+            let vec = assorted.iter_mut().find(|vec| vec[0].1.include(p))?;
+            vec.push((boundary, poly_boundary));
+            Some(())
+        })?;
+    let res = assorted
+        .into_iter()
+        .map(|vec| vec.into_iter().map(|(x, _)| x).collect())
+        .collect();
+    Some(res)
+}
+
+// --- connect_open_boundaries ---
+
+fn connect_open_boundaries<C>(
+    boundaries: &mut Vec<Wire>,
+    mut new_boundaries: Vec<Wire>,
+    divisor: (usize, usize),
+    edges: &[Edge<C>],
+    mut open: Vec<Wire>,
+) -> Option<()> {
+    let (v0, v1) = divisor;
+    let take_front = closure_take_front(edges);
+    let find_index = move |wire: &Wire, v: usize| {
+        wire.iter()
+            .position(|edge_index| take_front(*edge_index) == v)
+    };
+    let (i, index) = open
+        .iter()
+        .enumerate()
+        .find_map(|(i, wire)| find_index(wire, v0).map(|index| (i, index)))?;
+    if index != 0 {
+        open.swap(0, 1);
+    }
+    let (mut open_v1, mut open_v0) = (open.pop()?, open.pop()?);
+    let j = find_index(&open_v1, v1)?;
+    open_v0.rotate_left(i);
+    open_v1.rotate_left(j);
+    let (k, added_wire) = new_boundaries
+        .iter_mut()
+        .find_map(|wire| find_index(wire, v0).map(|index| (index, wire)))?;
+    added_wire.rotate_left(k);
+    open_v0.append(added_wire);
+    let l = find_index(&open_v0, v1)?;
+    open_v0.rotate_left(l);
+    open_v0.append(&mut open_v1);
+    *added_wire = open_v0;
+    *boundaries = new_boundaries;
+    Some(())
+}
+
+// --- common functions ---
+
+fn closure_take_front<C>(edges: &[Edge<C>]) -> impl Fn(EdgeIndex) -> usize + '_ {
+    move |edge_index: EdgeIndex| {
+        let index = edge_index.index;
+        let (v0, v1) = edges[index].vertices;
+        match edge_index.orientation {
+            true => v0,
+            false => v1,
+        }
+    }
+}
+
+fn closure_take_back<C>(edges: &[Edge<C>]) -> impl Fn(EdgeIndex) -> usize + '_ {
+    move |edge_index: EdgeIndex| {
+        let (v0, v1) = edges[edge_index.index].vertices;
+        match edge_index.orientation {
+            true => v1,
+            false => v0,
+        }
+    }
+}
+
+fn closure_get_poly<P: Clone>(
+    poly_edges: &[PolylineCurve<P>],
+) -> impl Fn(EdgeIndex) -> PolylineCurve<P> + '_ {
+    move |edge_index: EdgeIndex| {
+        let poly = &poly_edges[edge_index.index];
+        match edge_index.orientation {
+            true => poly.clone(),
+            false => poly.inverse(),
+        }
+    }
+}
+
+fn closure_to_poly<P, C>(tol: f64) -> impl Fn(&Edge<C>) -> PolylineCurve<P> + 'static
+where C: BoundedCurve + ParameterDivision1D<Point = P> {
+    move |Edge { curve, .. }: &Edge<C>| PolylineCurve::from_curve(curve, curve.range_tuple(), tol)
+}
+
+type ZippedEdge<'a> = (&'a EdgeIndex, &'a PolylineCurve<Point2>);
+type ZippedWire<'a> = (&'a Wire, &'a Vec<PolylineCurve<Point2>>);
+
+fn zip_boundaries<'a>(
+    boundaries: &'a [Wire],
+    param_boundaries: &'a [Vec<PolylineCurve<Point2>>],
+) -> impl Iterator<Item = ZippedEdge<'a>> {
+    boundaries
+        .iter()
+        .zip(param_boundaries)
+        .flat_map(move |(b, pb)| b.iter().zip(pb))
+}
+
+#[derive(Clone, Copy, Debug)]
+enum BoundaryType {
+    Positive,
+    Negative,
+    NotClosed,
+}
+
+fn boundary_type(poly: &PolylineCurve<Point2>) -> BoundaryType {
+    if poly[0].distance2(poly[poly.len() - 1]) < 1.0e-3 {
+        match poly.area() > 0.0 {
+            true => BoundaryType::Positive,
+            false => BoundaryType::Negative,
+        }
+    } else {
+        BoundaryType::NotClosed
+    }
+}
+
+fn periodic_iterator(
+    p: Point2,
+    (up, vp): (Option<f64>, Option<f64>),
+) -> impl Iterator<Item = Point2> {
+    let up_range = match up {
+        Some(up) => vec![-2.0 * up, -up, 0.0, up, 2.0 * up],
+        None => vec![0.0],
+    };
+    let vp_range = match vp {
+        Some(vp) => vec![-2.0 * vp, -vp, 0.0, vp, 2.0 * vp],
+        None => vec![0.0],
+    };
+    itertools::iproduct!(up_range, vp_range).map(move |(dx, dy)| p + Vector2::new(dx, dy))
+}
+
+fn nearest_correction(
+    (v0, _v1): (usize, usize),
+    wire0: &[EdgeIndex],
+    wire1: &[EdgeIndex],
+    vertices: &[Point3],
+    take_front: impl Fn(EdgeIndex) -> usize,
+) -> Option<(usize, usize)> {
+    let p0 = vertices[v0];
+    let v1 = wire1.iter().copied().map(&take_front).min_by(|v, w| {
+        p0.distance2(vertices[*v])
+            .partial_cmp(&p0.distance2(vertices[*w]))
+            .unwrap()
+    })?;
+    let p1 = vertices[v1];
+    let v0 = wire0.iter().copied().map(&take_front).min_by(|v, w| {
+        p1.distance2(vertices[*v])
+            .partial_cmp(&p1.distance2(vertices[*w]))
+            .unwrap()
+    })?;
+    Some((v0, v1))
+}
@@ -0,0 +1,37 @@
+//! Crate for operationg shapes. Provides boolean operations to Solid, and shape healing for importing shapes from other CAD systems.
+
+// SolveCraft: warnings are not errors in the vendored copy (newer compilers flag more).
+#![allow(dead_code)]
+#![warn(clippy::all, rust_2018_idioms)]
+#![warn(
+    missing_docs,
+    missing_debug_implementations,
+    trivial_casts,
+    trivial_numeric_casts,
+    unsafe_code,
+    unstable_features,
+    unused_import_braces,
+    unused_qualifications
+)]
+
+/// SolveCraft: report where an operation gives up (set `SHAPEOPS_TRACE`).
+macro_rules! tr {
+    ($e:expr, $tag:expr) => {
+        match $e {
+            Some(v) => v,
+            None => {
+                if std::env::var_os("SHAPEOPS_TRACE").is_some() {
+                    eprintln!("shapeops: gave up at {} ({}:{})", $tag, file!(), line!());
+                }
+                return None;
+            }
+        }
+    };
+}
+pub(crate) use tr;
+
+mod healing;
+pub use healing::{RobustSplitClosedEdgesAndFaces, SplitClosedEdgesAndFaces};
+mod transversal;
+pub use transversal::{and, or, ShapeOpsCurve, ShapeOpsSurface};
+mod alternative;
@@ -0,0 +1,140 @@
+#![allow(clippy::many_single_char_names)]
+
+use super::faces_classification::FacesClassification;
+use super::loops_store::*;
+use rustc_hash::FxHashMap as HashMap;
+use std::ops::Deref;
+use truck_meshalgo::prelude::*;
+use truck_topology::*;
+
+fn create_parameter_boundary<P, C, S>(
+    face: &Face<P, C, S>,
+    wire: &Wire<P, C>,
+    polys: &mut HashMap<EdgeID<C>, PolylineCurve<P>>,
+    tol: f64,
+) -> Option<PolylineCurve<Point2>>
+where
+    P: Copy,
+    C: BoundedCurve<Point = P> + ParameterDivision1D<Point = P>,
+    S: Clone + SearchParameter<D2, Point = P> + SearchNearestParameter<D2, Point = P>,
+{
+    let surface = face.surface();
+    let pt = wire.front_vertex().unwrap().point();
+    // SolveCraft: intersection-curve points lie on the surface only within the tolerance; fall
+    // back to the nearest parameter when the exact search fails.
+    let p: Point2 = surface
+        .search_parameter(pt, None, 100)
+        .or_else(|| surface.search_nearest_parameter(pt, None, 100))?
+        .into();
+    let vec = wire.edge_iter().try_fold(vec![p], |mut vec, edge| {
+        let poly = polys.entry(edge.id()).or_insert_with(|| {
+            let curve = edge.curve();
+            let div = curve.parameter_division(curve.range_tuple(), tol).1;
+            PolylineCurve(div)
+        });
+        let mut p = *vec.last().unwrap();
+        let closure = |q: &P| -> Option<Point2> {
+            p = surface
+                .search_parameter(*q, Some(p.into()), 100)
+                .or_else(|| surface.search_nearest_parameter(*q, Some(p.into()), 100))?
+                .into();
+            Some(p)
+        };
+        let add: Option<Vec<Point2>> = match edge.orientation() {
+            true => poly.iter().skip(1).map(closure).collect(),
+            false => poly.iter().rev().skip(1).map(closure).collect(),
+        };
+        vec.append(&mut add?);
+        Some(vec)
+    })?;
+    Some(PolylineCurve(vec))
+}
+
+#[derive(Clone, Debug)]
+struct WireChunk<'a, C> {
+    poly: PolylineCurve<Point2>,
+    wire: &'a BoundaryWire<Point3, C>,
+}
+
+type FaceWithShapesOpStatus<C, S> = (Face<Point3, C, S>, ShapesOpStatus);
+fn divide_one_face<C, S>(
+    face: &Face<Point3, C, S>,
+    loops: &Loops<Point3, C>,
+    tol: f64,
+) -> Option<Vec<FaceWithShapesOpStatus<C, S>>>
+where
+    C: BoundedCurve<Point = Point3> + ParameterDivision1D<Point = Point3>,
+    S: Clone + SearchParameter<D2, Point = Point3> + SearchNearestParameter<D2, Point = Point3>,
+{
+    let (mut pre_faces, mut negative_wires) = (Vec::new(), Vec::new());
+    let mut map = HashMap::default();
+    loops.iter().try_for_each(|wire| {
+        let poly = crate::tr!(create_parameter_boundary(face, wire, &mut map, tol), "parameter boundary");
+        match poly.area() > 0.0 {
+            true => pre_faces.push(vec![WireChunk { poly, wire }]),
+            false => negative_wires.push(WireChunk { poly, wire }),
+        }
+        Some(())
+    })?;
+    negative_wires.into_iter().try_for_each(|chunk| {
+        let pt = chunk.poly.front();
+        let op = crate::tr!(pre_faces.iter_mut().find(|face| face[0].poly.include(pt)), "loop not inside any face piece");
+        op.push(chunk);
+        Some(())
+    })?;
+    let vec: Vec<_> = pre_faces
+        .into_iter()
+        .map(|pre_face| {
+            let surface = face.surface();
+            let op = pre_face
+                .iter()
+                .find(|chunk| chunk.wire.status() != ShapesOpStatus::Unknown);
+            let status = match op {
+                Some(chunk) => chunk.wire.status(),
+                None => ShapesOpStatus::Unknown,
+            };
+            let wires: Vec<Wire<Point3, C>> = pre_face
+                .into_iter()
+                .map(|chunk| chunk.wire.deref().clone())
+                .collect();
+            let mut new_face = Face::debug_new(wires, surface);
+            if !face.orientation() {
+                new_face.invert();
+            }
+            (new_face, status)
+        })
+        .collect();
+    Some(vec)
+}
+
+pub fn divide_faces<C, S>(
+    shell: &Shell<Point3, C, S>,
+    loops_store: &LoopsStore<Point3, C>,
+    tol: f64,
+) -> Option<FacesClassification<Point3, C, S>>
+where
+    C: BoundedCurve<Point = Point3> + ParameterDivision1D<Point = Point3>,
+    S: Clone + SearchParameter<D2, Point = Point3> + SearchNearestParameter<D2, Point = Point3>,
+{
+    let mut res = FacesClassification::<Point3, C, S>::default();
+    shell
+        .iter()
+        .zip(loops_store)
+        .try_for_each(|(face, loops)| {
+            if loops
+                .iter()
+                .all(|wire| wire.status() == ShapesOpStatus::Unknown)
+            {
+                res.push(face.clone(), ShapesOpStatus::Unknown);
+            } else {
+                let vec = crate::tr!(divide_one_face(face, loops, tol), "divide one face");
+                vec.into_iter()
+                    .for_each(|(face, status)| res.push(face, status));
+            }
+            Some(())
+        })?;
+    Some(res)
+}
+
+#[cfg(test)]
+mod tests;
@@ -0,0 +1,226 @@
+use super::*;
+use shell::ShellCondition;
+use truck_geometry::prelude::*;
+use truck_topology::Vertex;
+const TOL: f64 = 0.05;
+
+fn line(v0: &Vertex<Point3>, v1: &Vertex<Point3>) -> Edge<Point3, BSplineCurve<Point3>> {
+    let curve = BSplineCurve::new(KnotVec::bezier_knot(1), vec![v0.point(), v1.point()]);
+    Edge::new(v0, v1, curve)
+}
+
+fn parabola(
+    v0: &Vertex<Point3>,
+    v1: &Vertex<Point3>,
+    pt: Point3,
+) -> Edge<Point3, BSplineCurve<Point3>> {
+    let curve = BSplineCurve::new(KnotVec::bezier_knot(2), vec![v0.point(), pt, v1.point()]);
+    Edge::new(v0, v1, curve)
+}
+
+#[test]
+fn divide_plane_test() {
+    let v = Vertex::news([
+        Point3::new(0.0, 0.0, 0.0),
+        Point3::new(0.0, 4.0, 0.0),
+        Point3::new(-1.0, 1.0, 0.0),
+        Point3::new(-1.0, 3.0, 0.0),
+        Point3::new(1.0, 1.0, 0.0),
+        Point3::new(1.0, 3.0, 0.0),
+    ]);
+    let edge = vec![
+        parabola(&v[0], &v[1], Point3::new(-4.0, 2.0, 0.0)),
+        parabola(&v[0], &v[1], Point3::new(4.0, 2.0, 0.0)),
+        line(&v[0], &v[1]),
+        parabola(&v[2], &v[3], Point3::new(-3.0, 2.0, 0.0)),
+        parabola(&v[2], &v[3], Point3::new(-1.0, 2.0, 0.0)),
+        parabola(&v[4], &v[5], Point3::new(1.0, 2.0, 0.0)),
+        parabola(&v[4], &v[5], Point3::new(3.0, 2.0, 0.0)),
+    ];
+    let wire: Vec<Wire<_, _>> = vec![
+        vec![edge[1].clone(), edge[0].inverse()].into(),
+        vec![edge[2].clone(), edge[0].inverse()].into(),
+        vec![edge[1].clone(), edge[2].inverse()].into(),
+        vec![edge[3].clone(), edge[4].inverse()].into(),
+        vec![edge[5].clone(), edge[6].inverse()].into(),
+    ];
+    let face = Face::new(
+        vec![wire[0].clone(), wire[3].clone(), wire[4].clone()],
+        Plane::new(
+            Point3::origin(),
+            Point3::new(1.0, 0.0, 0.0),
+            Point3::new(0.0, 1.0, 0.0),
+        ),
+    );
+    let loops: Loops<_, _> = vec![
+        BoundaryWire::new(wire[1].clone(), ShapesOpStatus::Or),
+        BoundaryWire::new(wire[2].clone(), ShapesOpStatus::And),
+        BoundaryWire::new(wire[3].clone(), ShapesOpStatus::Unknown),
+        BoundaryWire::new(wire[4].clone(), ShapesOpStatus::Unknown),
+    ]
+    .into_iter()
+    .collect();
+    let res = divide_one_face(&face, &loops, 0.01).unwrap();
+    assert_eq!(res.len(), 2);
+    let (mut or, mut and) = (true, true);
+    for (face, status) in res {
+        let bdd = face.absolute_boundaries();
+        match status {
+            ShapesOpStatus::Or => {
+                assert_eq!(bdd.len(), 2);
+                assert!(bdd[0] == wire[1] || bdd[0] == wire[3]);
+                assert!(bdd[1] == wire[1] || bdd[1] == wire[3]);
+                assert_ne!(bdd[0], bdd[1]);
+                assert!(or);
+                or = false;
+            }
+            ShapesOpStatus::And => {
+                assert_eq!(bdd.len(), 2);
+                assert!(bdd[0] == wire[2] || bdd[0] == wire[4]);
+                assert!(bdd[1] == wire[2] || bdd[1] == wire[4]);
+                assert_ne!(bdd[0], bdd[1]);
+                assert!(and);
+                and = false;
+            }
+            _ => panic!("There must be no unknown!"),
+        }
+    }
+}
+
+type AlternativeIntersection = crate::alternative::Alternative<
+    NurbsCurve<Vector4>,
+    IntersectionCurve<PolylineCurve<Point3>, AlternativeSurface>,
+>;
+type AlternativeSurface = crate::alternative::Alternative<BSplineSurface<Point3>, Plane>;
+
+crate::impl_from!(
+    NurbsCurve<Vector4>,
+    IntersectionCurve<PolylineCurve<Point3>, AlternativeSurface>
+);
+crate::impl_from!(BSplineSurface<Point3>, Plane);
+
+fn parabola_surfaces() -> (AlternativeSurface, AlternativeSurface) {
+    // define surfaces
+    #[rustfmt::skip]
+	let ctrl0 = vec![
+		vec![Point3::new(-1.0, -1.0, 3.0), Point3::new(-1.0, 0.0, -1.0), Point3::new(-1.0, 1.0, 3.0)],
+		vec![Point3::new(0.0, -1.0, -1.0), Point3::new(0.0, 0.0, -5.0), Point3::new(0.0, 1.0, -1.0)],
+		vec![Point3::new(1.0, -1.0, 3.0), Point3::new(1.0, 0.0, -1.0), Point3::new(1.0, 1.0, 3.0)],
+	];
+    #[rustfmt::skip]
+	let ctrl1 = vec![
+		vec![Point3::new(-1.0, -1.0, -3.0), Point3::new(-1.0, 0.0, 1.0), Point3::new(-1.0, 1.0, -3.0)],
+		vec![Point3::new(0.0, -1.0, 1.0), Point3::new(0.0, 0.0, 5.0), Point3::new(0.0, 1.0, 1.0)],
+		vec![Point3::new(1.0, -1.0, -3.0), Point3::new(1.0, 0.0, 1.0), Point3::new(1.0, 1.0, -3.0)],
+	];
+    (
+        BSplineSurface::new((KnotVec::bezier_knot(2), KnotVec::bezier_knot(2)), ctrl0).into(),
+        BSplineSurface::new((KnotVec::bezier_knot(2), KnotVec::bezier_knot(2)), ctrl1).into(),
+    )
+}
+
+#[test]
+fn independent_intersection() {
+    // prepare geoetries
+    let arc00: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new(
+        KnotVec::bezier_knot(2),
+        vec![
+            Vector4::new(1.0, 0.0, 1.0, 1.0),
+            Vector4::new(0.0, 1.0, 0.0, 0.0),
+            Vector4::new(-1.0, 0.0, 1.0, 1.0),
+        ],
+    ))
+    .into();
+    let arc01: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new(
+        KnotVec::bezier_knot(2),
+        vec![
+            Vector4::new(-1.0, 0.0, 1.0, 1.0),
+            Vector4::new(0.0, -1.0, 0.0, 0.0),
+            Vector4::new(1.0, 0.0, 1.0, 1.0),
+        ],
+    ))
+    .into();
+    let arc10: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new(
+        KnotVec::bezier_knot(2),
+        vec![
+            Vector4::new(1.0, 0.0, -1.0, 1.0),
+            Vector4::new(0.0, 1.0, 0.0, 0.0),
+            Vector4::new(-1.0, 0.0, -1.0, 1.0),
+        ],
+    ))
+    .into();
+    let arc11: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new(
+        KnotVec::bezier_knot(2),
+        vec![
+            Vector4::new(-1.0, 0.0, -1.0, 1.0),
+            Vector4::new(0.0, -1.0, 0.0, 0.0),
+            Vector4::new(1.0, 0.0, -1.0, 1.0),
+        ],
+    ))
+    .into();
+    let (surface0, surface1) = parabola_surfaces();
+    let plane0: AlternativeSurface = Plane::new(
+        Point3::new(0.0, 0.0, 1.0),
+        Point3::new(1.0, 0.0, 1.0),
+        Point3::new(0.0, 1.0, 1.0),
+    )
+    .into();
+    let plane1: AlternativeSurface = Plane::new(
+        Point3::new(0.0, 0.0, -1.0),
+        Point3::new(1.0, 0.0, -1.0),
+        Point3::new(0.0, 1.0, -1.0),
+    )
+    .into();
+
+    // prepare topologies
+    let v00 = Vertex::new(Point3::new(1.0, 0.0, 1.0));
+    let v01 = Vertex::new(Point3::new(-1.0, 0.0, 1.0));
+    let v10 = Vertex::new(Point3::new(1.0, 0.0, -1.0));
+    let v11 = Vertex::new(Point3::new(-1.0, 0.0, -1.0));
+    let wire0: Wire<_, _> = vec![Edge::new(&v00, &v01, arc00), Edge::new(&v01, &v00, arc01)].into();
+    let wire1: Wire<_, _> = vec![Edge::new(&v10, &v11, arc10), Edge::new(&v11, &v10, arc11)].into();
+    let shell0: Shell<_, _, _> = vec![
+        Face::new(vec![wire0.clone()], plane0),
+        Face::new(vec![wire0], surface0).inverse(),
+    ]
+    .into();
+    assert_eq!(shell0.shell_condition(), ShellCondition::Closed);
+    let shell1: Shell<_, _, _> = vec![
+        Face::new(vec![wire1.clone()], plane1).inverse(),
+        Face::new(vec![wire1], surface1),
+    ]
+    .into();
+    assert_eq!(shell1.shell_condition(), ShellCondition::Closed);
+    let poly_shell0 = shell0.triangulation(TOL);
+    let poly_shell1 = shell1.triangulation(TOL);
+
+    let LoopsStoreQuadruple {
+        geom_loops_store0: loops_store0,
+        geom_loops_store1: loops_store1,
+        ..
+    } = create_loops_stores(&shell0, &poly_shell0, &shell1, &poly_shell1, TOL).unwrap();
+    let [and0, or0, unknown0] = divide_faces(&shell0, &loops_store0, TOL)
+        .unwrap()
+        .and_or_unknown();
+    let [and1, or1, unknown1] = divide_faces(&shell1, &loops_store1, TOL)
+        .unwrap()
+        .and_or_unknown();
+    assert_eq!(and0.len(), 1);
+    assert_eq!(or0.len(), 1);
+    assert_eq!(unknown0.len(), 1);
+    assert_eq!(and1.len(), 1);
+    assert_eq!(or1.len(), 1);
+    assert_eq!(unknown1.len(), 1);
+
+    match unknown0[0].surface() {
+        AlternativeSurface::FirstType(_) => panic!("This face must plane!"),
+        AlternativeSurface::SecondType(_) => {}
+    }
+    match unknown1[0].surface() {
+        AlternativeSurface::FirstType(_) => panic!("This face must plane!"),
+        AlternativeSurface::SecondType(_) => {}
+    }
+
+    let and_shell: Shell<_, _, _> = vec![and0[0].clone(), and1[0].clone()].into();
+    assert_eq!(and_shell.shell_condition(), ShellCondition::Closed);
+}
@@ -0,0 +1,67 @@
+use super::loops_store::ShapesOpStatus;
+use rustc_hash::FxHashMap as HashMap;
+use truck_topology::*;
+
+#[derive(Clone, Debug)]
+pub struct FacesClassification<P, C, S> {
+    shell: Shell<P, C, S>,
+    status: HashMap<FaceID<S>, ShapesOpStatus>,
+}
+
+impl<P, C, S> Default for FacesClassification<P, C, S> {
+    fn default() -> Self {
+        Self {
+            shell: Default::default(),
+            status: HashMap::default(),
+        }
+    }
+}
+
+impl<P, C, S> FacesClassification<P, C, S> {
+    pub fn push(&mut self, face: Face<P, C, S>, status: ShapesOpStatus) {
+        self.status.insert(face.id(), status);
+        self.shell.push(face);
+    }
+
+    pub fn and_or_unknown(&self) -> [Shell<P, C, S>; 3] {
+        let [mut and, mut or, mut unknown] = <[Shell<P, C, S>; 3]>::default();
+        for face in &self.shell {
+            match self.status.get(&face.id()).unwrap() {
+                ShapesOpStatus::And => and.push(face.clone()),
+                ShapesOpStatus::Or => or.push(face.clone()),
+                ShapesOpStatus::Unknown => unknown.push(face.clone()),
+            }
+        }
+        [and, or, unknown]
+    }
+
+    pub fn integrate_by_component(&mut self) {
+        let [and, or, unknown] = self.and_or_unknown();
+        let and_boundary = and.extract_boundaries();
+        let or_boundary = or.extract_boundaries();
+        let components = unknown.connected_components();
+        for comp in components {
+            let boundary = comp.extract_boundaries();
+            if and_boundary
+                .iter()
+                .flatten()
+                .any(|edge| edge.id() == boundary[0][0].id())
+            {
+                comp.iter().for_each(|face| {
+                    *self.status.get_mut(&face.id()).unwrap() = ShapesOpStatus::And;
+                })
+            } else if or_boundary
+                .iter()
+                .flatten()
+                .any(|edge| edge.id() == boundary[0][0].id())
+            {
+                comp.iter().for_each(|face| {
+                    *self.status.get_mut(&face.id()).unwrap() = ShapesOpStatus::Or;
+                })
+            }
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests;
@@ -0,0 +1,131 @@
+use super::{
+    super::{divide_face, loops_store},
+    *,
+};
+use truck_geometry::prelude::*;
+use truck_meshalgo::prelude::*;
+use truck_topology::{shell::ShellCondition, Vertex};
+const TOL: f64 = 0.05;
+
+type AlternativeIntersection = crate::alternative::Alternative<
+    NurbsCurve<Vector4>,
+    IntersectionCurve<PolylineCurve<Point3>, AlternativeSurface>,
+>;
+type AlternativeSurface = crate::alternative::Alternative<BSplineSurface<Point3>, Plane>;
+
+fn parabola_surfaces() -> (AlternativeSurface, AlternativeSurface) {
+    // define surfaces
+    #[rustfmt::skip]
+	let ctrl0 = vec![
+		vec![Point3::new(-1.0, -1.0, 3.0), Point3::new(-1.0, 0.0, -1.0), Point3::new(-1.0, 1.0, 3.0)],
+		vec![Point3::new(0.0, -1.0, -1.0), Point3::new(0.0, 0.0, -5.0), Point3::new(0.0, 1.0, -1.0)],
+		vec![Point3::new(1.0, -1.0, 3.0), Point3::new(1.0, 0.0, -1.0), Point3::new(1.0, 1.0, 3.0)],
+	];
+    #[rustfmt::skip]
+	let ctrl1 = vec![
+		vec![Point3::new(-1.0, -1.0, -3.0), Point3::new(-1.0, 0.0, 1.0), Point3::new(-1.0, 1.0, -3.0)],
+		vec![Point3::new(0.0, -1.0, 1.0), Point3::new(0.0, 0.0, 5.0), Point3::new(0.0, 1.0, 1.0)],
+		vec![Point3::new(1.0, -1.0, -3.0), Point3::new(1.0, 0.0, 1.0), Point3::new(1.0, 1.0, -3.0)],
+	];
+    (
+        BSplineSurface::new((KnotVec::bezier_knot(2), KnotVec::bezier_knot(2)), ctrl0).into(),
+        BSplineSurface::new((KnotVec::bezier_knot(2), KnotVec::bezier_knot(2)), ctrl1).into(),
+    )
+}
+
+#[test]
+fn independent_intersection() {
+    // prepare geoetries
+    let arc00: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new(
+        KnotVec::bezier_knot(2),
+        vec![
+            Vector4::new(1.0, 0.0, 1.0, 1.0),
+            Vector4::new(0.0, 1.0, 0.0, 0.0),
+            Vector4::new(-1.0, 0.0, 1.0, 1.0),
+        ],
+    ))
+    .into();
+    let arc01: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new(
+        KnotVec::bezier_knot(2),
+        vec![
+            Vector4::new(-1.0, 0.0, 1.0, 1.0),
+            Vector4::new(0.0, -1.0, 0.0, 0.0),
+            Vector4::new(1.0, 0.0, 1.0, 1.0),
+        ],
+    ))
+    .into();
+    let arc10: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new(
+        KnotVec::bezier_knot(2),
+        vec![
+            Vector4::new(1.0, 0.0, -1.0, 1.0),
+            Vector4::new(0.0, 1.0, 0.0, 0.0),
+            Vector4::new(-1.0, 0.0, -1.0, 1.0),
+        ],
+    ))
+    .into();
+    let arc11: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new(
+        KnotVec::bezier_knot(2),
+        vec![
+            Vector4::new(-1.0, 0.0, -1.0, 1.0),
+            Vector4::new(0.0, -1.0, 0.0, 0.0),
+            Vector4::new(1.0, 0.0, -1.0, 1.0),
+        ],
+    ))
+    .into();
+    let (surface0, surface1) = parabola_surfaces();
+    let plane0: AlternativeSurface = Plane::new(
+        Point3::new(0.0, 0.0, 1.0),
+        Point3::new(1.0, 0.0, 1.0),
+        Point3::new(0.0, 1.0, 1.0),
+    )
+    .into();
+    let plane1: AlternativeSurface = Plane::new(
+        Point3::new(0.0, 0.0, -1.0),
+        Point3::new(1.0, 0.0, -1.0),
+        Point3::new(0.0, 1.0, -1.0),
+    )
+    .into();
+
+    // prepare topologies
+    let v00 = Vertex::new(Point3::new(1.0, 0.0, 1.0));
+    let v01 = Vertex::new(Point3::new(-1.0, 0.0, 1.0));
+    let v10 = Vertex::new(Point3::new(1.0, 0.0, -1.0));
+    let v11 = Vertex::new(Point3::new(-1.0, 0.0, -1.0));
+    let wire0: Wire<_, _> = vec![Edge::new(&v00, &v01, arc00), Edge::new(&v01, &v00, arc01)].into();
+    let wire1: Wire<_, _> = vec![Edge::new(&v10, &v11, arc10), Edge::new(&v11, &v10, arc11)].into();
+    let shell0: Shell<_, _, _> = vec![
+        Face::new(vec![wire0.clone()], plane0),
+        Face::new(vec![wire0], surface0).inverse(),
+    ]
+    .into();
+    assert_eq!(shell0.shell_condition(), ShellCondition::Closed);
+    let shell1: Shell<_, _, _> = vec![
+        Face::new(vec![wire1.clone()], plane1).inverse(),
+        Face::new(vec![wire1], surface1),
+    ]
+    .into();
+    assert_eq!(shell1.shell_condition(), ShellCondition::Closed);
+    let poly_shell0 = shell0.triangulation(TOL);
+    let poly_shell1 = shell1.triangulation(TOL);
+
+    let loops_store::LoopsStoreQuadruple {
+        geom_loops_store0: loops_store0,
+        geom_loops_store1: loops_store1,
+        ..
+    } = loops_store::create_loops_stores(&shell0, &poly_shell0, &shell1, &poly_shell1, TOL)
+        .unwrap();
+    let mut cls0 = divide_face::divide_faces(&shell0, &loops_store0, TOL).unwrap();
+    cls0.integrate_by_component();
+    let mut cls1 = divide_face::divide_faces(&shell1, &loops_store1, TOL).unwrap();
+    cls1.integrate_by_component();
+
+    let [mut and, mut or, _] = cls0.and_or_unknown();
+    let [and1, or1, _] = cls1.and_or_unknown();
+    and.extend(and1);
+    or.extend(or1);
+
+    assert_eq!(and.len(), 2);
+    assert_eq!(or.len(), 4);
+    assert_eq!(and.shell_condition(), ShellCondition::Closed);
+    assert_eq!(or.shell_condition(), ShellCondition::Closed);
+}
@@ -0,0 +1,151 @@
+use super::*;
+use truck_geometry::prelude::IntersectionCurve;
+use truck_meshalgo::prelude::*;
+use truck_topology::*;
+
+/// Only solids consisting of faces whose surface is implemented this trait can be used for set operations.
+pub trait ShapeOpsSurface:
+    ParametricSurface3D
+    + ParameterDivision2D
+    + SearchParameter<D2, Point = Point3>
+    + SearchNearestParameter<D2, Point = Point3>
+    + Invertible
+    + Send
+    + Sync {
+}
+impl<S> ShapeOpsSurface for S where S: ParametricSurface3D
+        + ParameterDivision2D
+        + SearchParameter<D2, Point = Point3>
+        + SearchNearestParameter<D2, Point = Point3>
+        + Invertible
+        + Send
+        + Sync
+{
+}
+
+/// Only solids consisting of edges whose curve is implemented this trait can be used for set operations.
+pub trait ShapeOpsCurve<S: ShapeOpsSurface>:
+    ParametricCurve3D
+    + ParameterDivision1D<Point = Point3>
+    + Cut
+    + Invertible
+    + From<IntersectionCurve<PolylineCurve<Point3>, S>>
+    + SearchParameter<D1, Point = Point3>
+    + SearchNearestParameter<D1, Point = Point3>
+    + Send
+    + Sync {
+}
+impl<C, S: ShapeOpsSurface> ShapeOpsCurve<S> for C where C: ParametricCurve3D
+        + ParameterDivision1D<Point = Point3>
+        + Cut
+        + Invertible
+        + From<IntersectionCurve<PolylineCurve<Point3>, S>>
+        + SearchParameter<D1, Point = Point3>
+        + SearchNearestParameter<D1, Point = Point3>
+        + Send
+        + Sync
+{
+}
+
+fn process_one_pair_of_shells<C: ShapeOpsCurve<S>, S: ShapeOpsSurface>(
+    shell0: &Shell<Point3, C, S>,
+    shell1: &Shell<Point3, C, S>,
+    tol: f64,
+) -> Option<[Shell<Point3, C, S>; 2]> {
+    nonpositive_tolerance!(tol);
+    // SolveCraft: robust meshing (boundaries from earlier operations lie on surfaces only
+    // within tolerance).
+    let poly_shell0 = shell0.robust_triangulation(tol);
+    let poly_shell1 = shell1.robust_triangulation(tol);
+    let loops_store::LoopsStoreQuadruple {
+        geom_loops_store0: loops_store0,
+        geom_loops_store1: loops_store1,
+        ..
+    } = crate::tr!(loops_store::create_loops_stores(shell0, &poly_shell0, shell1, &poly_shell1, tol), "loops stores");
+    let mut cls0 = crate::tr!(divide_face::divide_faces(shell0, &loops_store0, tol), "divide faces 0");
+    cls0.integrate_by_component();
+    let mut cls1 = crate::tr!(divide_face::divide_faces(shell1, &loops_store1, tol), "divide faces 1");
+    cls1.integrate_by_component();
+    let [mut and0, mut or0, unknown0] = cls0.and_or_unknown();
+    unknown0.into_iter().try_for_each(|face| {
+        let pt = face.boundaries()[0].vertex_iter().next().unwrap().point();
+        let dir = hash::take_one_unit(pt);
+        let count = poly_shell1.iter().try_fold(0, |count, face| {
+            let poly = face.surface()?;
+            Some(count + poly.signed_crossing_faces(pt, dir))
+        })?;
+        if count >= 1 {
+            and0.push(face);
+        } else {
+            or0.push(face);
+        }
+        Some(())
+    })?;
+    let [mut and1, mut or1, unknown1] = cls1.and_or_unknown();
+    unknown1.into_iter().try_for_each(|face| {
+        let pt = face.boundaries()[0].vertex_iter().next().unwrap().point();
+        let dir = hash::take_one_unit(pt);
+        let count = poly_shell0.iter().try_fold(0, |count, face| {
+            let poly = face.surface()?;
+            Some(count + poly.signed_crossing_faces(pt, dir))
+        })?;
+        if count >= 1 {
+            and1.push(face);
+        } else {
+            or1.push(face);
+        }
+        Some(())
+    })?;
+    and0.append(&mut and1);
+    or0.append(&mut or1);
+    Some([and0, or0])
+}
+
+/// AND operation between two solids.
+pub fn and<C: ShapeOpsCurve<S>, S: ShapeOpsSurface>(
+    solid0: &Solid<Point3, C, S>,
+    solid1: &Solid<Point3, C, S>,
+    tol: f64,
+) -> Option<Solid<Point3, C, S>> {
+    let mut iter0 = solid0.boundaries().iter();
+    let mut iter1 = solid1.boundaries().iter();
+    let shell0 = iter0.next().unwrap();
+    let shell1 = iter1.next().unwrap();
+    let [mut and_shell, _] = process_one_pair_of_shells(shell0, shell1, tol)?;
+    for shell in iter0 {
+        let [res, _] = process_one_pair_of_shells(&and_shell, shell, tol)?;
+        and_shell = res;
+    }
+    for shell in iter1 {
+        let [res, _] = process_one_pair_of_shells(&and_shell, shell, tol)?;
+        and_shell = res;
+    }
+    let boundaries = and_shell.connected_components();
+    Some(Solid::new(boundaries))
+}
+
+/// OR operation between two solids.
+pub fn or<C: ShapeOpsCurve<S>, S: ShapeOpsSurface>(
+    solid0: &Solid<Point3, C, S>,
+    solid1: &Solid<Point3, C, S>,
+    tol: f64,
+) -> Option<Solid<Point3, C, S>> {
+    let mut iter0 = solid0.boundaries().iter();
+    let mut iter1 = solid1.boundaries().iter();
+    let shell0 = iter0.next().unwrap();
+    let shell1 = iter1.next().unwrap();
+    let [_, mut or_shell] = process_one_pair_of_shells(shell0, shell1, tol)?;
+    for shell in iter0 {
+        let [_, res] = process_one_pair_of_shells(&or_shell, shell, tol)?;
+        or_shell = res;
+    }
+    for shell in iter1 {
+        let [_, res] = process_one_pair_of_shells(&or_shell, shell, tol)?;
+        or_shell = res;
+    }
+    let boundaries = or_shell.connected_components();
+    Some(Solid::new(boundaries))
+}
+
+#[cfg(test)]
+mod tests;
@@ -0,0 +1,21 @@
+use truck_meshalgo::prelude::*;
+use truck_modeling::*;
+
+#[test]
+fn punched_cube() {
+    let v = builder::vertex(Point3::origin());
+    let e = builder::tsweep(&v, Vector3::unit_x());
+    let f = builder::tsweep(&e, Vector3::unit_y());
+    let cube = builder::tsweep(&f, Vector3::unit_z());
+
+    let v = builder::vertex(Point3::new(0.5, 0.25, -0.5));
+    let w = builder::rsweep(&v, Point3::new(0.5, 0.5, 0.0), Vector3::unit_z(), Rad(7.0));
+    let f = builder::try_attach_plane(&[w]).unwrap();
+    let mut cylinder = builder::tsweep(&f, Vector3::unit_z() * 2.0);
+    cylinder.not();
+    let and = crate::and(&cube, &cylinder, 0.05).unwrap();
+
+    let poly = and.triangulation(0.01).to_polygon();
+    let file = std::fs::File::create("punched-cube.obj").unwrap();
+    obj::write(&poly, file).unwrap();
+}
@@ -0,0 +1,174 @@
+use truck_base::cgmath64::*;
+use truck_geometry::prelude::*;
+use truck_meshalgo::prelude::*;
+
+/// polyline base ntersection curve with parameter
+#[derive(Debug, Clone, derive_more::Deref, derive_more::DerefMut)]
+pub struct IntersectionCurveWithParameters<S> {
+    #[deref]
+    #[deref_mut]
+    ic: IntersectionCurve<PolylineCurve<Point3>, S>,
+    params0: PolylineCurve<Point2>,
+    params1: PolylineCurve<Point2>,
+}
+
+impl<S> From<IntersectionCurveWithParameters<S>> for IntersectionCurve<PolylineCurve<Point3>, S> {
+    fn from(a: IntersectionCurveWithParameters<S>) -> Self { a.ic }
+}
+
+impl<S> IntersectionCurveWithParameters<S>
+where S: ParametricSurface3D + SearchNearestParameter<D2, Point = Point3>
+{
+    pub fn try_new(
+        surface0: S,
+        surface1: S,
+        poly: PolylineCurve<Point3>,
+        tol: f64,
+    ) -> Option<Self> {
+        let mut polyline = PolylineCurve(Vec::new());
+        let mut params0 = PolylineCurve(Vec::new());
+        let mut params1 = PolylineCurve(Vec::new());
+        for p in poly.windows(2) {
+            let n = (p[1] - p[0]).normalize();
+            let (q, p0, p1) = double_projection(&surface0, None, &surface1, None, p[0], n, 100)?;
+            polyline.push(q);
+            params0.push(p0);
+            params1.push(p1);
+        }
+        let (q, p0, p1) = if poly[0].near(&poly[poly.len() - 1]) {
+            (polyline[0], params0[0], params1[0])
+        } else {
+            let n = (poly[poly.len() - 1] - poly[poly.len() - 2]).normalize();
+            double_projection(
+                &surface0,
+                None,
+                &surface1,
+                None,
+                poly[poly.len() - 1],
+                n,
+                100,
+            )?
+        };
+        polyline.push(q);
+        params0.push(p0);
+        params1.push(p1);
+        Some(Self {
+            ic: IntersectionCurve::new_unchecked(
+                Box::new(surface0),
+                Box::new(surface1),
+                polyline,
+                tol,
+            ),
+            params0,
+            params1,
+        })
+    }
+}
+
+impl<S> ParametricCurve for IntersectionCurveWithParameters<S>
+where S: ParametricSurface3D + SearchNearestParameter<D2, Point = Point3>
+{
+    type Point = Point3;
+    type Vector = Vector3;
+    #[inline(always)]
+    fn subs(&self, t: f64) -> Point3 { self.ic.subs(t) }
+    #[inline(always)]
+    fn der(&self, t: f64) -> Vector3 { self.ic.der(t) }
+    #[inline(always)]
+    fn der2(&self, t: f64) -> Vector3 { self.ic.der2(t) }
+    #[inline(always)]
+    fn parameter_range(&self) -> ParameterRange { self.ic.parameter_range() }
+}
+
+impl<S> BoundedCurve for IntersectionCurveWithParameters<S> where S: ParametricSurface3D + SearchNearestParameter<D2, Point = Point3> {}
+
+impl<S> ParameterDivision1D for IntersectionCurveWithParameters<S>
+where S: ParametricSurface3D + SearchNearestParameter<D2, Point = Point3>
+{
+    type Point = Point3;
+    #[inline(always)]
+    fn parameter_division(&self, range: (f64, f64), tol: f64) -> (Vec<f64>, Vec<Point3>) {
+        self.ic.parameter_division(range, tol)
+    }
+}
+
+impl<S> Cut for IntersectionCurveWithParameters<S>
+where S: ParametricSurface3D + SearchNearestParameter<D2, Point = Point3>
+{
+    #[inline(always)]
+    fn cut(&mut self, t: f64) -> Self {
+        Self {
+            ic: self.ic.cut(t),
+            params0: self.params0.cut(t),
+            params1: self.params1.cut(t),
+        }
+    }
+}
+
+impl<S: Clone> Invertible for IntersectionCurveWithParameters<S> {
+    fn invert(&mut self) {
+        self.ic.invert();
+        self.params0.invert();
+        self.params1.invert();
+    }
+}
+
+impl<S> SearchParameter<D1> for IntersectionCurveWithParameters<S>
+where S: ParametricSurface3D + SearchNearestParameter<D2, Point = Point3>
+{
+    type Point = Point3;
+    fn search_parameter<H: Into<SPHint1D>>(
+        &self,
+        point: Point3,
+        hint: H,
+        trials: usize,
+    ) -> Option<f64> {
+        self.ic.search_parameter(point, hint, trials)
+    }
+}
+
+impl<S> SearchNearestParameter<D1> for IntersectionCurveWithParameters<S>
+where S: ParametricSurface3D + SearchNearestParameter<D2, Point = Point3>
+{
+    type Point = Point3;
+    fn search_nearest_parameter<H: Into<SPHint1D>>(
+        &self,
+        point: Point3,
+        hint: H,
+        trials: usize,
+    ) -> Option<f64> {
+        self.ic.search_nearest_parameter(point, hint, trials)
+    }
+}
+
+pub fn intersection_curves<S>(
+    surface0: S,
+    polygon0: &PolygonMesh,
+    surface1: S,
+    polygon1: &PolygonMesh,
+    tol: f64,
+) -> Vec<(
+    PolylineCurve<Point3>,
+    Option<IntersectionCurveWithParameters<S>>,
+)>
+where
+    S: ParametricSurface3D + SearchNearestParameter<D2, Point = Point3>,
+{
+    let interferences = polygon0.extract_interference(polygon1);
+    let polylines = super::polyline_construction::construct_polylines(&interferences);
+    polylines
+        .into_iter()
+        .map(|polyline| {
+            let curve = IntersectionCurveWithParameters::try_new(
+                surface0.clone(),
+                surface1.clone(),
+                polyline.clone(),
+                tol,
+            );
+            (polyline, curve)
+        })
+        .collect()
+}
+
+#[cfg(test)]
+mod tests;
@@ -0,0 +1,173 @@
+use super::*;
+use std::f64::consts::PI;
+
+mod double_projection_tests {
+    use super::*;
+
+    fn create_axis(n: Vector3) -> (Vector3, Vector3) {
+        let idx = if n[0].abs() < n[1].abs() { 0 } else { 1 };
+        let idx = if n[idx].abs() < n[2].abs() { idx } else { 2 };
+        let mut e = Vector3::zero();
+        e[idx] = 1.0;
+        let x = n.cross(e).normalize();
+        (x, n.cross(x))
+    }
+
+    fn exec_plane_case() {
+        let c0 = Point3::new(
+            2.0 * rand::random::<f64>() - 1.0,
+            2.0 * rand::random::<f64>() - 1.0,
+            2.0 * rand::random::<f64>() - 1.0,
+        );
+        let n0 = Vector3::new(
+            2.0 * rand::random::<f64>() - 1.0,
+            2.0 * rand::random::<f64>() - 1.0,
+            2.0 * rand::random::<f64>() - 1.0,
+        )
+        .normalize();
+        let (x, y) = create_axis(n0);
+        let plane0 = Plane::new(c0, c0 + x, c0 + y);
+        let c1 = Point3::new(
+            2.0 * rand::random::<f64>() - 1.0,
+            2.0 * rand::random::<f64>() - 1.0,
+            2.0 * rand::random::<f64>() - 1.0,
+        );
+        let n1 = Vector3::new(
+            2.0 * rand::random::<f64>() - 1.0,
+            2.0 * rand::random::<f64>() - 1.0,
+            2.0 * rand::random::<f64>() - 1.0,
+        )
+        .normalize();
+        let (x, y) = create_axis(n1);
+        let plane1 = Plane::new(c1, c1 + x, c1 + y);
+        let n = n0.cross(n1).normalize();
+        let mut o = None;
+        for i in 0..10 {
+            let t = i as f64;
+            let p = Point3::origin() + t * n;
+            let (q, p0, p1) = double_projection(&plane0, None, &plane1, None, p, n, 100)
+                .unwrap_or_else(|| panic!("plane0: {plane0:?}\nplane1: {plane1:?}\n p: {p:?}"));
+            assert_near!(q, plane0.subs(p0.x, p0.y));
+            assert_near!(q, plane1.subs(p1.x, p1.y));
+            if let Some(o) = o {
+                assert_near!(q.distance2(o), t * t);
+            } else {
+                o = Some(q);
+            }
+        }
+    }
+
+    #[test]
+    fn plane_case() { (0..100).for_each(|_| exec_plane_case()); }
+
+    #[test]
+    fn sphere_case() {
+        let sphere0 = Sphere::new(Point3::new(0.0, 0.0, 1.0), f64::sqrt(2.0));
+        let sphere1 = Sphere::new(Point3::new(0.0, 0.0, -1.0), f64::sqrt(2.0));
+        for _ in 0..100 {
+            let t = 2.0 * PI * rand::random::<f64>();
+            let r = 2.0 * rand::random::<f64>();
+            let p = Point3::new(r * f64::cos(t), r * f64::sin(t), 0.0);
+            let n = Vector3::new(-f64::sin(t), f64::cos(t), 0.0);
+            let (q, p0, p1) = double_projection(&sphere0, None, &sphere1, None, p, n, 100)
+                .unwrap_or_else(|| panic!("p: {p:?}"));
+            assert_near!(q, sphere0.subs(p0.x, p0.y));
+            assert_near!(q, sphere1.subs(p1.x, p1.y));
+            assert_near!(q, Point3::new(f64::cos(t), f64::sin(t), 0.0));
+        }
+    }
+}
+
+#[test]
+fn intersection_curve_sphere_case() {
+    let sphere0 = Sphere::new(Point3::new(0.0, 0.0, 1.0), f64::sqrt(2.0));
+    let sphere1 = Sphere::new(Point3::new(0.0, 0.0, -1.0), f64::sqrt(2.0));
+    const M: usize = 5;
+    let polyline = (0..=M)
+        .map(|i| {
+            let t = 2.0 * PI * i as f64 / M as f64;
+            Point3::new(0.8 * f64::cos(t), 0.8 * f64::sin(t), 0.0)
+        })
+        .collect::<PolylineCurve<_>>();
+    let curve: IntersectionCurve<_, _> =
+        IntersectionCurveWithParameters::try_new(sphere0, sphere1, polyline, 0.5)
+            .unwrap()
+            .into();
+
+    const N: usize = 100;
+    let mut sum = 0.0;
+    let (t0, t1) = curve.range_tuple();
+    for i in 0..N {
+        let t = t0 + (t1 - t0) * i as f64 / N as f64;
+        let pt = curve.subs(t);
+        assert_near!(pt.distance(Point3::origin()), 1.0);
+        let vec = curve.der(t);
+        assert!(pt.dot(vec).so_small(), "{i} {t} {vec:?}");
+        assert!(vec[2].so_small());
+        sum += vec.magnitude() * (t1 - t0) / N as f64;
+    }
+    assert!(f64::abs(sum - 2.0 * PI) < 0.1, "{}", sum);
+
+    let theta = 2.0 * PI * rand::random::<f64>();
+    let pt = Point3::new(f64::cos(theta), f64::sin(theta), 0.0);
+    let t = curve.search_parameter(pt, None, 1).unwrap();
+    assert!(curve.subs(t).near(&pt));
+    let pt = Point3::new(1.1 * f64::cos(theta), 1.1 * f64::sin(theta), 0.0);
+    assert!(curve.search_parameter(pt, None, 1).is_none());
+
+    let mut curve0 = curve.clone();
+    let curve1 = curve0.cut(2.5);
+    assert_near!(curve0.front(), curve.front());
+    assert_near!(curve0.back(), curve.subs(2.5));
+    assert_near!(curve1.front(), curve.subs(2.5));
+    assert_near!(curve1.back(), curve.back());
+    let mut curve0 = curve.clone();
+    let curve1 = curve0.cut(2.0);
+    assert_near!(curve0.front(), curve.front());
+    assert_near!(curve0.back(), curve.subs(2.0));
+    assert_near!(curve1.front(), curve.subs(2.0));
+    assert_near!(curve1.back(), curve.back());
+}
+
+#[test]
+fn collide_parabola() {
+    const TOL: f64 = 0.05;
+
+    // define surfaces
+    #[rustfmt::skip]
+	let ctrl0 = vec![
+		vec![Point3::new(-1.0, -1.0, 3.0), Point3::new(-1.0, 0.0, -1.0), Point3::new(-1.0, 1.0, 3.0)],
+		vec![Point3::new(0.0, -1.0, -1.0), Point3::new(0.0, 0.0, -5.0), Point3::new(0.0, 1.0, -1.0)],
+		vec![Point3::new(1.0, -1.0, 3.0), Point3::new(1.0, 0.0, -1.0), Point3::new(1.0, 1.0, 3.0)],
+	];
+    #[rustfmt::skip]
+	let ctrl1 = vec![
+		vec![Point3::new(-1.0, -1.0, -3.0), Point3::new(-1.0, 0.0, 1.0), Point3::new(-1.0, 1.0, -3.0)],
+		vec![Point3::new(0.0, -1.0, 1.0), Point3::new(0.0, 0.0, 5.0), Point3::new(0.0, 1.0, 1.0)],
+		vec![Point3::new(1.0, -1.0, -3.0), Point3::new(1.0, 0.0, 1.0), Point3::new(1.0, 1.0, -3.0)],
+	];
+    let surface0 = BSplineSurface::new((KnotVec::bezier_knot(2), KnotVec::bezier_knot(2)), ctrl0);
+    let surface1 = BSplineSurface::new((KnotVec::bezier_knot(2), KnotVec::bezier_knot(2)), ctrl1);
+
+    // meshing surface
+    let instant = std::time::Instant::now();
+    let polygon0 = StructuredMesh::from_surface(&surface0, surface0.range_tuple(), TOL).destruct();
+    let polygon1 = StructuredMesh::from_surface(&surface1, surface1.range_tuple(), TOL).destruct();
+    println!("Meshing Surfaces: {}s", instant.elapsed().as_secs_f64());
+    // extract intersection curves
+    let instant = std::time::Instant::now();
+    let curves = intersection_curves(surface0, &polygon0, surface1, &polygon1, TOL);
+    println!(
+        "Extracting Intersection: {}s",
+        instant.elapsed().as_secs_f64()
+    );
+    assert_eq!(curves.len(), 1);
+    let curve = curves[0].1.clone().unwrap();
+    const N: usize = 100;
+    for i in 0..N {
+        let t1 = curve.range_tuple().1;
+        let t = t1 * i as f64 / N as f64;
+        let pt = curve.subs(t);
+        assert_near!(pt.distance(Point3::origin()) * 0.5, f64::sqrt(0.5) * 0.5);
+    }
+}
@@ -0,0 +1,579 @@
+#![allow(clippy::many_single_char_names)]
+
+use super::*;
+use rustc_hash::FxHashMap as HashMap;
+use truck_base::cgmath64::*;
+use truck_geometry::prelude::IntersectionCurve;
+use truck_meshalgo::prelude::*;
+use truck_topology::{Vertex, *};
+
+type PolylineCurve = truck_meshalgo::prelude::PolylineCurve<Point3>;
+
+#[derive(Clone, Copy, Debug, PartialEq, Eq)]
+pub enum ShapesOpStatus {
+    Unknown,
+    And,
+    Or,
+}
+
+impl ShapesOpStatus {
+    fn not(self) -> Self {
+        match self {
+            Self::Unknown => Self::Unknown,
+            Self::And => Self::Or,
+            Self::Or => Self::And,
+        }
+    }
+}
+
+#[derive(Clone, Debug)]
+pub struct BoundaryWire<P, C> {
+    wire: Wire<P, C>,
+    status: ShapesOpStatus,
+}
+
+impl<P, C> BoundaryWire<P, C> {
+    #[inline(always)]
+    pub fn new(wire: Wire<P, C>, status: ShapesOpStatus) -> Self { Self { wire, status } }
+    #[inline(always)]
+    pub fn status(&self) -> ShapesOpStatus { self.status }
+    #[inline(always)]
+    pub fn invert(&mut self) {
+        self.wire.invert();
+        self.status = self.status.not();
+    }
+    #[inline(always)]
+    pub fn inverse(&self) -> Self {
+        Self {
+            wire: self.wire.inverse(),
+            status: self.status.not(),
+        }
+    }
+}
+
+impl ShapesOpStatus {
+    fn from_is_curve<C, S>(curve: &IntersectionCurve<C, S>) -> Option<ShapesOpStatus>
+    where
+        C: ParametricCurve3D + BoundedCurve,
+        S: ParametricSurface3D + SearchNearestParameter<D2, Point = Point3>, {
+        let (t0, t1) = curve.range_tuple();
+        let t = (t0 + t1) / 2.0;
+        let (_, pt0, pt1) = curve.search_triple(t)?;
+        let der = curve.leader().der(t);
+        let normal0 = curve.surface0().normal(pt0[0], pt0[1]);
+        let normal1 = curve.surface1().normal(pt1[0], pt1[1]);
+        match normal0.cross(der).dot(normal1) > 0.0 {
+            true => Some(ShapesOpStatus::Or),
+            false => Some(ShapesOpStatus::And),
+        }
+    }
+}
+
+impl<P, C> std::ops::Deref for BoundaryWire<P, C> {
+    type Target = Wire<P, C>;
+    #[inline(always)]
+    fn deref(&self) -> &Self::Target { &self.wire }
+}
+
+impl<P, C> std::ops::DerefMut for BoundaryWire<P, C> {
+    #[inline(always)]
+    fn deref_mut(&mut self) -> &mut Self::Target { &mut self.wire }
+}
+
+#[derive(Clone, Debug)]
+pub struct Loops<P, C>(Vec<BoundaryWire<P, C>>);
+#[derive(Clone, Debug)]
+pub struct LoopsStore<P, C>(Vec<Loops<P, C>>);
+
+impl<P, C> std::ops::Deref for Loops<P, C> {
+    type Target = Vec<BoundaryWire<P, C>>;
+    #[inline(always)]
+    fn deref(&self) -> &Self::Target { &self.0 }
+}
+
+impl<P, C> std::ops::DerefMut for Loops<P, C> {
+    #[inline(always)]
+    fn deref_mut(&mut self) -> &mut Self::Target { &mut self.0 }
+}
+
+impl<P, C> std::ops::Deref for LoopsStore<P, C> {
+    type Target = Vec<Loops<P, C>>;
+    #[inline(always)]
+    fn deref(&self) -> &Self::Target { &self.0 }
+}
+
+impl<P, C> std::ops::DerefMut for LoopsStore<P, C> {
+    #[inline(always)]
+    fn deref_mut(&mut self) -> &mut Self::Target { &mut self.0 }
+}
+
+impl<P, C> FromIterator<BoundaryWire<P, C>> for Loops<P, C> {
+    #[inline(always)]
+    fn from_iter<I: IntoIterator<Item = BoundaryWire<P, C>>>(iter: I) -> Self {
+        Self(Vec::from_iter(iter))
+    }
+}
+
+impl<'a, P, C, S> From<&'a Face<P, C, S>> for Loops<P, C> {
+    #[inline(always)]
+    fn from(face: &'a Face<P, C, S>) -> Loops<P, C> {
+        face.absolute_boundaries()
+            .iter()
+            .map(|wire| BoundaryWire::new(wire.clone(), ShapesOpStatus::Unknown))
+            .collect()
+    }
+}
+
+impl<'a, P: 'a, C: 'a, S: 'a> FromIterator<&'a Face<P, C, S>> for LoopsStore<P, C> {
+    fn from_iter<I: IntoIterator<Item = &'a Face<P, C, S>>>(iter: I) -> Self {
+        Self(iter.into_iter().map(Loops::from).collect())
+    }
+}
+
+impl<'a, P, C> IntoIterator for &'a LoopsStore<P, C> {
+    type Item = <&'a Vec<Loops<P, C>> as IntoIterator>::Item;
+    type IntoIter = <&'a Vec<Loops<P, C>> as IntoIterator>::IntoIter;
+    fn into_iter(self) -> Self::IntoIter { self.0.iter() }
+}
+
+#[derive(Clone, Debug, Copy, PartialEq)]
+enum ParameterKind {
+    Front,
+    Back,
+    Inner(f64),
+}
+
+impl ParameterKind {
+    fn try_new(t: f64, (t0, t1): (f64, f64)) -> Option<ParameterKind> {
+        if t0.near(&t) {
+            Some(ParameterKind::Front)
+        } else if t1.near(&t) {
+            Some(ParameterKind::Back)
+        } else if t0 < t && t < t1 {
+            Some(ParameterKind::Inner(t))
+        } else {
+            None
+        }
+    }
+}
+
+impl<P: Copy, C: Clone> Loops<P, C> {
+    fn search_parameter(&self, pt: P) -> Option<(usize, usize, ParameterKind)>
+    where C: BoundedCurve<Point = P> + SearchParameter<D1, Point = P> {
+        self.iter()
+            .enumerate()
+            .flat_map(move |(i, wire)| wire.iter().enumerate().map(move |(j, edge)| (i, j, edge)))
+            .find_map(|(i, j, edge)| {
+                let curve = edge.curve();
+                curve.search_parameter(pt, None, 1).and_then(|t| {
+                    let kind = ParameterKind::try_new(t, curve.range_tuple())?;
+                    Some((i, j, kind))
+                })
+            })
+    }
+
+    fn change_vertex(
+        &mut self,
+        old_vertex: &Vertex<P>,
+        new_vertex: &Vertex<P>,
+        emap: &mut HashMap<EdgeID<C>, Edge<P, C>>,
+    ) {
+        self.iter_mut()
+            .flat_map(|wire| wire.iter_mut())
+            .for_each(|edge| {
+                let mut new_edge = if edge.absolute_front() == old_vertex {
+                    emap.entry(edge.id()).or_insert_with(|| {
+                        Edge::new(new_vertex, edge.absolute_back(), edge.curve())
+                    })
+                } else if edge.absolute_back() == old_vertex {
+                    emap.entry(edge.id()).or_insert_with(|| {
+                        Edge::new(edge.absolute_front(), new_vertex, edge.curve())
+                    })
+                } else {
+                    return;
+                }
+                .clone();
+                if !edge.orientation() {
+                    new_edge.invert();
+                }
+                // Remove the edge from the HashMap when it is no longer there because ID reassignment will occur.
+                if edge.count() == 1 {
+                    emap.remove(&edge.id());
+                }
+                *edge = new_edge;
+            })
+    }
+
+    fn swap_edge_into_wire(&mut self, edge_id: EdgeID<C>, new_wire: &Wire<P, C>) {
+        self.iter_mut().for_each(|wire| {
+            let mut iter = wire.iter().enumerate();
+            if let Some((idx, edge)) = iter.find(|(_, edge)| edge.id() == edge_id) {
+                if edge.orientation() {
+                    wire.swap_edge_into_wire(idx, new_wire.clone());
+                } else {
+                    wire.swap_edge_into_wire(idx, new_wire.inverse());
+                }
+            }
+        });
+    }
+
+    #[inline(always)]
+    fn add_independent_loop(&mut self, r#loop: BoundaryWire<P, C>) {
+        self.push(r#loop.inverse());
+        self.push(r#loop);
+    }
+
+    fn add_edge(
+        &mut self,
+        edge0: Edge<P, C>,
+        status: ShapesOpStatus,
+    ) -> [Option<(usize, usize)>; 2] {
+        let a = self.iter().enumerate().find_map(|(i, wire)| {
+            wire.iter().enumerate().find_map(|(j, edge)| {
+                if edge.front() == edge0.back() {
+                    Some((i, j))
+                } else {
+                    None
+                }
+            })
+        });
+        let b = self.iter().enumerate().find_map(|(i, wire)| {
+            wire.iter().enumerate().find_map(|(j, edge)| {
+                if edge.front() == edge0.front() {
+                    Some((i, j))
+                } else {
+                    None
+                }
+            })
+        });
+        if let Some((wire_index0, edge_index0)) = a {
+            self[wire_index0].rotate_left(edge_index0);
+            self[wire_index0].push_front(edge0.clone());
+            self[wire_index0].push_back(edge0.inverse());
+        }
+        match (a, b) {
+            (Some((wire_index0, edge_index0)), Some((wire_index1, edge_index1))) => {
+                if wire_index0 == wire_index1 {
+                    let len = self[wire_index0].len() - 2;
+                    let edge_index1 = (len + edge_index1 - edge_index0) % len + 1;
+                    let new_wire = self[wire_index0].split_off(edge_index1);
+                    self[wire_index0].status = status;
+                    self.push(BoundaryWire::new(new_wire, status.not()));
+                } else {
+                    let mut new_wire0 = self[wire_index1].clone();
+                    let mut new_wire1 = new_wire0.split_off(edge_index1);
+                    new_wire0.append(&mut self[wire_index0]);
+                    new_wire0.append(&mut new_wire1);
+                    self[wire_index0] = new_wire0;
+                    self.swap_remove(wire_index1);
+                }
+            }
+            (None, Some((wire_index1, edge_index1))) => {
+                self[wire_index1].rotate_left(edge_index1);
+                self[wire_index1].push_front(edge0.inverse());
+                self[wire_index1].push_back(edge0);
+            }
+            (None, None) => self.push(BoundaryWire::new(
+                vec![edge0.inverse(), edge0].into(),
+                ShapesOpStatus::Unknown,
+            )),
+            _ => {}
+        }
+        [a, b]
+    }
+}
+
+impl<P: Copy + Tolerance, C: Clone> LoopsStore<P, C> {
+    #[inline(always)]
+    fn change_vertex(
+        &mut self,
+        old_vertex: &Vertex<P>,
+        new_vertex: &Vertex<P>,
+        emap: &mut HashMap<EdgeID<C>, Edge<P, C>>,
+    ) {
+        self.iter_mut()
+            .for_each(|loops| loops.change_vertex(old_vertex, new_vertex, emap));
+    }
+
+    #[inline(always)]
+    fn swap_edge_into_wire(&mut self, edge_id: EdgeID<C>, new_wire: &Wire<P, C>) {
+        self.iter_mut()
+            .for_each(|loops| loops.swap_edge_into_wire(edge_id, new_wire))
+    }
+
+    fn add_polygon_vertex(
+        &mut self,
+        loops_index: usize,
+        v: &Vertex<P>,
+        emap: &mut HashMap<EdgeID<C>, Edge<P, C>>,
+    ) -> Option<(usize, usize, ParameterKind)>
+    where
+        C: Cut<Point = P> + SearchParameter<D1, Point = P>,
+    {
+        let pt = v.point();
+        let (wire_index, edge_index, kind) = self[loops_index].search_parameter(pt)?;
+        match kind {
+            ParameterKind::Front => {
+                let old_vertex = self[loops_index][wire_index][edge_index]
+                    .absolute_front()
+                    .clone();
+                self.change_vertex(&old_vertex, v, emap);
+            }
+            ParameterKind::Back => {
+                let old_vertex = self[loops_index][wire_index][edge_index]
+                    .absolute_back()
+                    .clone();
+                self.change_vertex(&old_vertex, v, emap);
+            }
+            ParameterKind::Inner(t) => {
+                let edge = self[loops_index][wire_index][edge_index].absolute_clone();
+                let edge_id = edge.id();
+                let (edge0, edge1) = edge.cut_with_parameter(v, t)?;
+                let new_wire: Wire<_, _> = vec![edge0, edge1].into();
+                self.swap_edge_into_wire(edge_id, &new_wire);
+            }
+        }
+        Some((wire_index, edge_index, kind))
+    }
+}
+
+impl<C> LoopsStore<Point3, C> {
+    fn add_geom_vertex<S>(
+        &mut self,
+        (loops_index, wire_index, edge_index): (usize, usize, usize),
+        v: &Vertex<Point3>,
+        kind: ParameterKind,
+        another_surface: &S,
+        emap: &mut HashMap<EdgeID<C>, Edge<Point3, C>>,
+    ) -> Option<()>
+    where
+        C: Cut<Point = Point3, Vector = Vector3> + SearchNearestParameter<D1, Point = Point3>,
+        S: ParametricSurface3D + SearchNearestParameter<D2, Point = Point3>,
+    {
+        match kind {
+            ParameterKind::Front => {
+                let old_vertex = self[loops_index][wire_index][edge_index]
+                    .absolute_front()
+                    .clone();
+                v.set_point(old_vertex.point());
+                self.change_vertex(&old_vertex, v, emap);
+            }
+            ParameterKind::Back => {
+                let old_vertex = self[loops_index][wire_index][edge_index]
+                    .absolute_back()
+                    .clone();
+                v.set_point(old_vertex.point());
+                self.change_vertex(&old_vertex, v, emap);
+            }
+            ParameterKind::Inner(_) => {
+                let curve = self[loops_index][wire_index][edge_index].curve();
+                let (pt, t, _) =
+                    curve_surface_projection(&curve, None, another_surface, None, v.point(), 100)?;
+                v.set_point(pt);
+                let edge = self[loops_index][wire_index][edge_index].absolute_clone();
+                let edge_id = edge.id();
+                let (edge0, edge1) = edge.cut_with_parameter(v, t)?;
+                let new_wire: Wire<_, _> = vec![edge0, edge1].into();
+                self.swap_edge_into_wire(edge_id, &new_wire);
+            }
+        }
+        Some(())
+    }
+}
+
+fn curve_surface_projection<C, S>(
+    curve: &C,
+    curve_hint: Option<f64>,
+    surface: &S,
+    surface_hint: Option<(f64, f64)>,
+    point: Point3,
+    trials: usize,
+) -> Option<(Point3, f64, Point2)>
+where
+    C: ParametricCurve3D + SearchNearestParameter<D1, Point = Point3>,
+    S: ParametricSurface3D + SearchNearestParameter<D2, Point = Point3>,
+{
+    if trials == 0 {
+        return None;
+    }
+    let t = curve.search_nearest_parameter(point, curve_hint, 10)?;
+    let pt0 = curve.subs(t);
+    let (u, v) = surface.search_nearest_parameter(point, surface_hint, 10)?;
+    let pt1 = surface.subs(u, v);
+    if point.near(&pt0) && point.near(&pt1) && pt0.near(&pt1) {
+        Some((point, t, Point2::new(u, v)))
+    } else {
+        let l = curve.der(t);
+        let n = surface.normal(u, v);
+        let t0 = (pt1 - pt0).dot(n) / l.dot(n);
+        curve_surface_projection(
+            curve,
+            Some(t),
+            surface,
+            Some((u, v)),
+            pt0 + t0 * l,
+            trials - 1,
+        )
+    }
+}
+
+fn create_independent_loop<P, C, D>(mut poly_curve0: C) -> Wire<P, D>
+where
+    C: Cut<Point = P>,
+    D: From<C>, {
+    let (t0, t1) = poly_curve0.range_tuple();
+    let t = (t0 + t1) / 2.0;
+    let poly_curve1 = poly_curve0.cut(t);
+    let v0 = Vertex::new(poly_curve0.front());
+    let v1 = Vertex::new(poly_curve1.front());
+    let edge0 = Edge::new(&v0, &v1, poly_curve0.into());
+    let edge1 = Edge::new(&v1, &v0, poly_curve1.into());
+    vec![edge0, edge1].into()
+}
+
+#[allow(dead_code)]
+pub struct LoopsStoreQuadruple<C> {
+    pub geom_loops_store0: LoopsStore<Point3, C>,
+    pub poly_loops_store0: LoopsStore<Point3, PolylineCurve>,
+    pub geom_loops_store1: LoopsStore<Point3, C>,
+    pub poly_loops_store1: LoopsStore<Point3, PolylineCurve>,
+}
+
+pub fn create_loops_stores<C, S>(
+    geom_shell0: &Shell<Point3, C, S>,
+    poly_shell0: &Shell<Point3, PolylineCurve, Option<PolygonMesh>>,
+    geom_shell1: &Shell<Point3, C, S>,
+    poly_shell1: &Shell<Point3, PolylineCurve, Option<PolygonMesh>>,
+    tol: f64,
+) -> Option<LoopsStoreQuadruple<C>>
+where
+    C: SearchNearestParameter<D1, Point = Point3>
+        + SearchParameter<D1, Point = Point3>
+        + Cut<Point = Point3, Vector = Vector3>
+        + From<IntersectionCurve<PolylineCurve, S>>,
+    S: ParametricSurface3D + SearchNearestParameter<D2, Point = Point3>,
+{
+    let mut geom_loops_store0: LoopsStore<_, _> = geom_shell0.face_iter().collect();
+    let mut poly_loops_store0: LoopsStore<_, _> = poly_shell0.face_iter().collect();
+    let mut geom_loops_store1: LoopsStore<_, _> = geom_shell1.face_iter().collect();
+    let mut poly_loops_store1: LoopsStore<_, _> = poly_shell1.face_iter().collect();
+    let store0_len = geom_loops_store0.len();
+    let store1_len = geom_loops_store1.len();
+    (0..store0_len)
+        .flat_map(move |i| (0..store1_len).map(move |j| (i, j)))
+        .try_for_each(|(face_index0, face_index1)| {
+            let ori0 = geom_shell0[face_index0].orientation();
+            let ori1 = geom_shell1[face_index1].orientation();
+            let surface0 = geom_shell0[face_index0].surface();
+            let surface1 = geom_shell1[face_index1].surface();
+            let polygon0 = crate::tr!(poly_shell0[face_index0].surface(), "polygon0");
+            let polygon1 = crate::tr!(poly_shell1[face_index1].surface(), "polygon1");
+            intersection_curve::intersection_curves(
+                surface0.clone(),
+                &polygon0,
+                surface1.clone(),
+                &polygon1,
+                tol,
+            )
+            .into_iter()
+            .try_for_each(|(polyline, intersection_curve)| {
+                let mut intersection_curve = crate::tr!(intersection_curve, "intersection curve").into();
+                let status = crate::tr!(ShapesOpStatus::from_is_curve(&intersection_curve), "status");
+                let (status0, status1) = match (ori0, ori1) {
+                    (true, true) => (status, status.not()),
+                    (true, false) => (status.not(), status.not()),
+                    (false, true) => (status, status),
+                    (false, false) => (status.not(), status),
+                };
+                if polyline.front().near(&polyline.back()) {
+                    let poly_wire = create_independent_loop(polyline);
+                    poly_loops_store0[face_index0]
+                        .add_independent_loop(BoundaryWire::new(poly_wire.clone(), status0));
+                    poly_loops_store1[face_index1]
+                        .add_independent_loop(BoundaryWire::new(poly_wire, status1));
+                    let geom_wire = create_independent_loop(intersection_curve);
+                    geom_loops_store0[face_index0]
+                        .add_independent_loop(BoundaryWire::new(geom_wire.clone(), status0));
+                    geom_loops_store1[face_index1]
+                        .add_independent_loop(BoundaryWire::new(geom_wire, status1));
+                } else {
+                    let pv0 = Vertex::new(polyline.front());
+                    let pv1 = Vertex::new(polyline.back());
+                    let gv0 = Vertex::new(polyline.front());
+                    let gv1 = Vertex::new(polyline.back());
+                    let mut pemap0 = HashMap::default();
+                    let mut pemap1 = HashMap::default();
+                    let mut gemap0 = HashMap::default();
+                    let mut gemap1 = HashMap::default();
+                    let idx00 =
+                        poly_loops_store0.add_polygon_vertex(face_index0, &pv0, &mut pemap0);
+                    if let Some((wire_index, edge_index, kind)) = idx00 {
+                        crate::tr!(geom_loops_store0.add_geom_vertex(
+                            (face_index0, wire_index, edge_index),
+                            &gv0,
+                            kind,
+                            &surface1,
+                            &mut gemap0,
+                        ), "add_geom_vertex 1");
+                        let polyline = intersection_curve.editor().leader;
+                        *polyline.first_mut().unwrap() = gv0.point();
+                    }
+                    let idx01 =
+                        poly_loops_store0.add_polygon_vertex(face_index0, &pv1, &mut pemap1);
+                    if let Some((wire_index, edge_index, kind)) = idx01 {
+                        crate::tr!(geom_loops_store0.add_geom_vertex(
+                            (face_index0, wire_index, edge_index),
+                            &gv1,
+                            kind,
+                            &surface1,
+                            &mut gemap1,
+                        ), "add_geom_vertex 2");
+                        let polyline = intersection_curve.editor().leader;
+                        *polyline.last_mut().unwrap() = gv1.point();
+                    }
+                    let idx10 =
+                        poly_loops_store1.add_polygon_vertex(face_index1, &pv0, &mut pemap0);
+                    if let Some((wire_index, edge_index, kind)) = idx10 {
+                        crate::tr!(geom_loops_store1.add_geom_vertex(
+                            (face_index1, wire_index, edge_index),
+                            &gv0,
+                            kind,
+                            &surface0,
+                            &mut gemap0,
+                        ), "add_geom_vertex 3");
+                        let polyline = intersection_curve.editor().leader;
+                        *polyline.first_mut().unwrap() = gv0.point();
+                    }
+                    let idx11 =
+                        poly_loops_store1.add_polygon_vertex(face_index1, &pv1, &mut pemap1);
+                    if let Some((wire_index, edge_index, kind)) = idx11 {
+                        crate::tr!(geom_loops_store1.add_geom_vertex(
+                            (face_index1, wire_index, edge_index),
+                            &gv1,
+                            kind,
+                            &surface0,
+                            &mut gemap1,
+                        ), "add_geom_vertex 4");
+                        let polyline = intersection_curve.editor().leader;
+                        *polyline.last_mut().unwrap() = gv1.point();
+                    }
+                    let pedge = Edge::new(&pv0, &pv1, polyline);
+                    let gedge = Edge::new(&gv0, &gv1, intersection_curve.into());
+                    poly_loops_store0[face_index0].add_edge(pedge.clone(), status0);
+                    geom_loops_store0[face_index0].add_edge(gedge.clone(), status0);
+                    poly_loops_store1[face_index1].add_edge(pedge, status1);
+                    geom_loops_store1[face_index1].add_edge(gedge, status1);
+                }
+                Some(())
+            })
+        })?;
+    Some(LoopsStoreQuadruple {
+        geom_loops_store0,
+        poly_loops_store0,
+        geom_loops_store1,
+        poly_loops_store1,
+    })
+}
+
+#[cfg(test)]
+mod tests;
@@ -0,0 +1,776 @@
+use super::*;
+use std::fmt::Debug;
+use truck_geometry::prelude::*;
+const TOL: f64 = 0.01;
+
+crate::impl_from!(
+    NurbsCurve<Vector4>,
+    IntersectionCurve<PolylineCurve, BSplineSurface<Point3>>
+);
+type AlternativeIntersection = crate::alternative::Alternative<
+    NurbsCurve<Vector4>,
+    IntersectionCurve<PolylineCurve, BSplineSurface<Point3>>,
+>;
+
+struct DebugDisplay<'a, T, Format> {
+    entity: &'a T,
+    format: Format,
+}
+
+impl<'a, P: Debug, C: Debug> Debug for DebugDisplay<'a, Loops<P, C>, WireDisplayFormat> {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        f.debug_tuple("Loops")
+            .field(
+                &self
+                    .entity
+                    .0
+                    .iter()
+                    .map(|wire| wire.display(self.format))
+                    .collect::<Vec<_>>(),
+            )
+            .finish()
+    }
+}
+
+impl<'a, P: Debug, C: Debug> Debug for DebugDisplay<'a, LoopsStore<P, C>, WireDisplayFormat> {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        f.debug_list()
+            .entries(self.entity.0.iter().map(|loops| DebugDisplay {
+                entity: loops,
+                format: self.format,
+            }))
+            .finish()
+    }
+}
+
+impl<P: Debug, C: Debug> LoopsStore<P, C> {
+    fn display(&self, format: WireDisplayFormat) -> DebugDisplay<'_, Self, WireDisplayFormat> {
+        DebugDisplay {
+            entity: self,
+            format,
+        }
+    }
+}
+
+fn parabola_surfaces() -> (BSplineSurface<Point3>, BSplineSurface<Point3>) {
+    // define surfaces
+    #[rustfmt::skip]
+	let ctrl0 = vec![
+		vec![Point3::new(-1.0, -1.0, 3.0), Point3::new(-1.0, 0.0, -1.0), Point3::new(-1.0, 1.0, 3.0)],
+		vec![Point3::new(0.0, -1.0, -1.0), Point3::new(0.0, 0.0, -5.0), Point3::new(0.0, 1.0, -1.0)],
+		vec![Point3::new(1.0, -1.0, 3.0), Point3::new(1.0, 0.0, -1.0), Point3::new(1.0, 1.0, 3.0)],
+	];
+    #[rustfmt::skip]
+	let ctrl1 = vec![
+		vec![Point3::new(-1.0, -1.0, -3.0), Point3::new(-1.0, 0.0, 1.0), Point3::new(-1.0, 1.0, -3.0)],
+		vec![Point3::new(0.0, -1.0, 1.0), Point3::new(0.0, 0.0, 5.0), Point3::new(0.0, 1.0, 1.0)],
+		vec![Point3::new(1.0, -1.0, -3.0), Point3::new(1.0, 0.0, 1.0), Point3::new(1.0, 1.0, -3.0)],
+	];
+    (
+        BSplineSurface::new((KnotVec::bezier_knot(2), KnotVec::bezier_knot(2)), ctrl0),
+        BSplineSurface::new((KnotVec::bezier_knot(2), KnotVec::bezier_knot(2)), ctrl1),
+    )
+}
+
+#[test]
+fn independent_intersection() {
+    // prepare geoetries
+    let arc00: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new(
+        KnotVec::bezier_knot(2),
+        vec![
+            Vector4::new(1.0, 0.0, 1.0, 1.0),
+            Vector4::new(0.0, 1.0, 0.0, 0.0),
+            Vector4::new(-1.0, 0.0, 1.0, 1.0),
+        ],
+    ))
+    .into();
+    let arc01 = NurbsCurve::new(BSplineCurve::new(
+        KnotVec::bezier_knot(2),
+        vec![
+            Vector4::new(-1.0, 0.0, 1.0, 1.0),
+            Vector4::new(0.0, -1.0, 0.0, 0.0),
+            Vector4::new(1.0, 0.0, 1.0, 1.0),
+        ],
+    ))
+    .into();
+    let arc10: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new(
+        KnotVec::bezier_knot(2),
+        vec![
+            Vector4::new(1.0, 0.0, -1.0, 1.0),
+            Vector4::new(0.0, 1.0, 0.0, 0.0),
+            Vector4::new(-1.0, 0.0, -1.0, 1.0),
+        ],
+    ))
+    .into();
+    let arc11 = NurbsCurve::new(BSplineCurve::new(
+        KnotVec::bezier_knot(2),
+        vec![
+            Vector4::new(-1.0, 0.0, -1.0, 1.0),
+            Vector4::new(0.0, -1.0, 0.0, 0.0),
+            Vector4::new(1.0, 0.0, -1.0, 1.0),
+        ],
+    ))
+    .into();
+    let (surface0, surface1) = parabola_surfaces();
+
+    // prepare topologies
+    let v00 = Vertex::new(Point3::new(1.0, 0.0, 1.0));
+    let v01 = Vertex::new(Point3::new(-1.0, 0.0, 1.0));
+    let v10 = Vertex::new(Point3::new(1.0, 0.0, -1.0));
+    let v11 = Vertex::new(Point3::new(-1.0, 0.0, -1.0));
+    let wire0: Wire<_, _> = vec![Edge::new(&v00, &v01, arc00), Edge::new(&v01, &v00, arc01)].into();
+    let wire1: Wire<_, _> = vec![Edge::new(&v10, &v11, arc10), Edge::new(&v11, &v10, arc11)].into();
+    let geom_shell0: Shell<_, _, _> = vec![Face::new(vec![wire0], surface0).inverse()].into();
+    let geom_shell1: Shell<_, _, _> = vec![Face::new(vec![wire1], surface1)].into();
+    let poly_shell0 = geom_shell0.triangulation(TOL);
+    let poly_shell1 = geom_shell1.triangulation(TOL);
+
+    // exec create loops stores!
+    let LoopsStoreQuadruple {
+        geom_loops_store0,
+        geom_loops_store1,
+        ..
+    } = create_loops_stores(&geom_shell0, &poly_shell0, &geom_shell1, &poly_shell1, TOL).unwrap();
+
+    // check the topology
+    let vertex_format = VertexDisplayFormat::AsPoint;
+    let edge_id_format = EdgeDisplayFormat::VerticesTupleAndID { vertex_format };
+    let wire_id_format = WireDisplayFormat::EdgesListTuple {
+        edge_format: edge_id_format,
+    };
+    let edge_geom_format = EdgeDisplayFormat::VerticesTupleAndCurve { vertex_format };
+    let wire_geom_format = WireDisplayFormat::EdgesListTuple {
+        edge_format: edge_geom_format,
+    };
+    assert_eq!(geom_loops_store0.len(), 1);
+    assert_eq!(geom_loops_store0[0].len(), 3);
+    assert_eq!(geom_loops_store0[0][0].len(), 2);
+    assert_eq!(geom_loops_store0[0][1].len(), 2);
+    assert_eq!(geom_loops_store0[0][2].len(), 2);
+    assert!(
+        geom_loops_store0[0][0].is_closed(),
+        "{:?}",
+        geom_loops_store0.display(wire_id_format)
+    );
+    assert!(
+        geom_loops_store0[0][1].is_closed(),
+        "{:?}",
+        geom_loops_store0.display(wire_id_format)
+    );
+    assert!(
+        geom_loops_store0[0][2].is_closed(),
+        "{:?}",
+        geom_loops_store0.display(wire_id_format)
+    );
+    assert!(
+        geom_loops_store0[0][0].is_geometric_consistent(),
+        "{:?}",
+        geom_loops_store0.display(wire_geom_format)
+    );
+    assert!(
+        geom_loops_store0[0][1].is_geometric_consistent(),
+        "{:?}",
+        geom_loops_store0.display(wire_geom_format)
+    );
+    assert!(
+        geom_loops_store0[0][2].is_geometric_consistent(),
+        "{:?}",
+        geom_loops_store0.display(wire_geom_format)
+    );
+    assert_eq!(geom_loops_store1.len(), 1);
+    assert_eq!(geom_loops_store1[0].len(), 3);
+    assert_eq!(geom_loops_store1[0][0].len(), 2);
+    assert_eq!(geom_loops_store1[0][1].len(), 2);
+    assert_eq!(geom_loops_store1[0][2].len(), 2);
+    assert!(
+        geom_loops_store1[0][0].is_closed(),
+        "{:?}",
+        geom_loops_store1.display(wire_id_format)
+    );
+    assert!(
+        geom_loops_store1[0][1].is_closed(),
+        "{:?}",
+        geom_loops_store1.display(wire_id_format)
+    );
+    assert!(
+        geom_loops_store1[0][2].is_closed(),
+        "{:?}",
+        geom_loops_store1.display(wire_id_format)
+    );
+    assert!(
+        geom_loops_store1[0][0].is_geometric_consistent(),
+        "{:?}",
+        geom_loops_store1.display(wire_geom_format)
+    );
+    assert!(
+        geom_loops_store1[0][1].is_geometric_consistent(),
+        "{:?}",
+        geom_loops_store1.display(wire_geom_format)
+    );
+    assert!(
+        geom_loops_store1[0][2].is_geometric_consistent(),
+        "{:?}",
+        geom_loops_store1.display(wire_geom_format)
+    );
+
+    // check the boundary status
+    let mut flags = [true; 3];
+    for i in 0..3 {
+        let bw = &geom_loops_store0[0][i];
+        match bw.status() {
+            ShapesOpStatus::Unknown => {
+                let curve = bw[0].oriented_curve();
+                assert_near!(curve.subs(0.0)[2], 1.0);
+                assert!(flags[0]);
+                flags[0] = false;
+            }
+            ShapesOpStatus::And => {
+                let curve = bw[0].oriented_curve();
+                let pt = curve.subs(0.5) - Point3::origin();
+                let der = curve.der(0.5);
+                assert!(pt.cross(der)[2] > 0.0);
+                assert!(flags[1]);
+                flags[1] = false;
+            }
+            ShapesOpStatus::Or => {
+                let curve = bw[0].oriented_curve();
+                let pt = curve.subs(0.5) - Point3::origin();
+                let der = curve.der(0.5);
+                assert!(pt.cross(der)[2] < 0.0);
+                assert!(flags[2]);
+                flags[2] = false;
+            }
+        }
+    }
+    let mut flags = [true; 3];
+    for i in 0..3 {
+        let bw = &geom_loops_store1[0][i];
+        match bw.status() {
+            ShapesOpStatus::Unknown => {
+                let curve = bw[0].oriented_curve();
+                assert_near!(curve.subs(0.0)[2], -1.0);
+                assert!(flags[0]);
+                flags[0] = false;
+            }
+            ShapesOpStatus::Or => {
+                let curve = bw[0].oriented_curve();
+                let pt = curve.subs(0.5) - Point3::origin();
+                let der = curve.der(0.5);
+                assert!(pt.cross(der)[2] < 0.0);
+                assert!(flags[1]);
+                flags[1] = false;
+            }
+            ShapesOpStatus::And => {
+                let curve = bw[0].oriented_curve();
+                let pt = curve.subs(0.5) - Point3::origin();
+                let der = curve.der(0.5);
+                assert!(pt.cross(der)[2] > 0.0);
+                assert!(flags[2]);
+                flags[2] = false;
+            }
+        }
+    }
+}
+
+#[test]
+fn rotated_intersection() {
+    // prepare geoetries
+    let arc00: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new(
+        KnotVec::bezier_knot(2),
+        vec![
+            Vector4::new(1.0, 0.0, 1.0, 1.0),
+            Vector4::new(0.0, 1.0, 0.0, 0.0),
+            Vector4::new(-1.0, 0.0, 1.0, 1.0),
+        ],
+    ))
+    .into();
+    let arc01: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new(
+        KnotVec::bezier_knot(2),
+        vec![
+            Vector4::new(-1.0, 0.0, 1.0, 1.0),
+            Vector4::new(0.0, -1.0, 0.0, 0.0),
+            Vector4::new(1.0, 0.0, 1.0, 1.0),
+        ],
+    ))
+    .into();
+    let arc02: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new(
+        KnotVec::bezier_knot(2),
+        vec![
+            Vector4::new(1.0, 0.0, 1.0, 1.0),
+            Vector4::new(0.0, 0.0, -3.0, 1.0),
+            Vector4::new(-1.0, 0.0, 1.0, 1.0),
+        ],
+    ))
+    .into();
+    let arc10: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new(
+        KnotVec::bezier_knot(2),
+        vec![
+            Vector4::new(0.0, -1.0, -1.0, 1.0),
+            Vector4::new(1.0, 0.0, 0.0, 0.0),
+            Vector4::new(0.0, 1.0, -1.0, 1.0),
+        ],
+    ))
+    .into();
+    let arc11: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new(
+        KnotVec::bezier_knot(2),
+        vec![
+            Vector4::new(0.0, 1.0, -1.0, 1.0),
+            Vector4::new(-1.0, 0.0, 0.0, 0.0),
+            Vector4::new(0.0, -1.0, -1.0, 1.0),
+        ],
+    ))
+    .into();
+    let arc12: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new(
+        KnotVec::bezier_knot(2),
+        vec![
+            Vector4::new(0.0, -1.0, -1.0, 1.0),
+            Vector4::new(0.0, 0.0, 3.0, 1.0),
+            Vector4::new(0.0, 1.0, -1.0, 1.0),
+        ],
+    ))
+    .into();
+    let (surface0, surface1) = parabola_surfaces();
+
+    let v00 = Vertex::new(Point3::new(1.0, 0.0, 1.0));
+    let v01 = Vertex::new(Point3::new(-1.0, 0.0, 1.0));
+    let edge00 = Edge::new(&v00, &v01, arc00);
+    let edge01 = Edge::new(&v01, &v00, arc01);
+    let edge02 = Edge::new(&v00, &v01, arc02);
+    let wire00: Wire<_, _> = vec![edge00, edge02.inverse()].into();
+    let wire01: Wire<_, _> = vec![edge01, edge02].into();
+    let face00 = Face::new(vec![wire00], surface0.clone());
+    let face01 = Face::new(vec![wire01], surface0);
+    let geom_shell0: Shell<_, _, _> = vec![face00.inverse(), face01.inverse()].into();
+
+    let v10 = Vertex::new(Point3::new(0.0, -1.0, -1.0));
+    let v11 = Vertex::new(Point3::new(0.0, 1.0, -1.0));
+    let edge10 = Edge::new(&v10, &v11, arc10);
+    let edge11 = Edge::new(&v11, &v10, arc11);
+    let edge12 = Edge::new(&v10, &v11, arc12);
+    let wire10: Wire<_, _> = vec![edge10, edge12.inverse()].into();
+    let wire11: Wire<_, _> = vec![edge12, edge11].into();
+    let face10 = Face::new(vec![wire10], surface1.clone());
+    let face11 = Face::new(vec![wire11], surface1);
+    let geom_shell1: Shell<_, _, _> = vec![face10, face11].into();
+
+    let poly_shell0 = geom_shell0.triangulation(TOL);
+    let poly_shell1 = geom_shell1.triangulation(TOL);
+    let file = std::fs::File::create("polyshell0.obj").unwrap();
+    obj::write(&poly_shell0.to_polygon(), file).unwrap();
+    let file = std::fs::File::create("polyshell1.obj").unwrap();
+    obj::write(&poly_shell1.to_polygon(), file).unwrap();
+    let LoopsStoreQuadruple {
+        geom_loops_store0,
+        geom_loops_store1,
+        ..
+    } = create_loops_stores(&geom_shell0, &poly_shell0, &geom_shell1, &poly_shell1, TOL).unwrap();
+
+    let vertex_format = VertexDisplayFormat::AsPoint;
+    let edge_id_format = EdgeDisplayFormat::VerticesTupleAndID { vertex_format };
+    let wire_id_format = WireDisplayFormat::EdgesListTuple {
+        edge_format: edge_id_format,
+    };
+    let edge_geom_format = EdgeDisplayFormat::VerticesTupleAndCurve { vertex_format };
+    let wire_geom_format = WireDisplayFormat::EdgesListTuple {
+        edge_format: edge_geom_format,
+    };
+    assert_eq!(geom_loops_store0.len(), 2);
+    assert_eq!(geom_loops_store0[0].len(), 2);
+    assert!(
+        geom_loops_store0[0][0].is_closed(),
+        "{:?}",
+        geom_loops_store0.display(wire_id_format)
+    );
+    assert!(
+        geom_loops_store0[0][1].is_closed(),
+        "{:?}",
+        geom_loops_store0.display(wire_id_format)
+    );
+    assert!(
+        geom_loops_store0[0][0].is_geometric_consistent(),
+        "{:?}",
+        geom_loops_store0.display(wire_geom_format)
+    );
+    assert!(
+        geom_loops_store0[0][1].is_geometric_consistent(),
+        "{:?}",
+        geom_loops_store0.display(wire_geom_format)
+    );
+    let (a, b) = (geom_loops_store0[0][0].len(), geom_loops_store0[0][1].len());
+    let compatible0 = match geom_loops_store0[0][0].status() {
+        ShapesOpStatus::And => a == 3,
+        ShapesOpStatus::Or => a == 5,
+        _ => false,
+    };
+    let compatible1 = match geom_loops_store0[0][1].status() {
+        ShapesOpStatus::And => b == 3,
+        ShapesOpStatus::Or => b == 5,
+        _ => false,
+    };
+    assert!(
+        compatible0 && compatible1 && a * b == 15,
+        "{:?}",
+        geom_loops_store0.display(wire_id_format)
+    );
+    assert_eq!(geom_loops_store0[1].len(), 2);
+    assert!(
+        geom_loops_store0[1][0].is_closed(),
+        "{:?}",
+        geom_loops_store0.display(wire_id_format)
+    );
+    assert!(
+        geom_loops_store0[1][1].is_closed(),
+        "{:?}",
+        geom_loops_store0.display(wire_id_format)
+    );
+    assert!(
+        geom_loops_store0[1][0].is_geometric_consistent(),
+        "{:?}",
+        geom_loops_store0.display(wire_geom_format)
+    );
+    assert!(
+        geom_loops_store0[1][1].is_geometric_consistent(),
+        "{:?}",
+        geom_loops_store0.display(wire_geom_format)
+    );
+    let (a, b) = (geom_loops_store0[1][0].len(), geom_loops_store0[1][1].len());
+    let compatible0 = match geom_loops_store0[1][0].status() {
+        ShapesOpStatus::And => a == 3,
+        ShapesOpStatus::Or => a == 5,
+        _ => false,
+    };
+    let compatible1 = match geom_loops_store0[1][1].status() {
+        ShapesOpStatus::And => b == 3,
+        ShapesOpStatus::Or => b == 5,
+        _ => false,
+    };
+    assert!(
+        compatible0 && compatible1 && a * b == 15,
+        "{:?}",
+        geom_loops_store0.display(wire_id_format)
+    );
+    assert_eq!(geom_loops_store1.len(), 2);
+    assert_eq!(geom_loops_store1[0].len(), 2);
+    assert!(
+        geom_loops_store1[0][0].is_closed(),
+        "{:?}",
+        geom_loops_store1.display(wire_id_format)
+    );
+    assert!(
+        geom_loops_store1[0][1].is_closed(),
+        "{:?}",
+        geom_loops_store1.display(wire_id_format)
+    );
+    assert!(
+        geom_loops_store1[0][0].is_geometric_consistent(),
+        "{:?}",
+        geom_loops_store1.display(wire_geom_format)
+    );
+    assert!(
+        geom_loops_store1[0][1].is_geometric_consistent(),
+        "{:?}",
+        geom_loops_store1.display(wire_geom_format)
+    );
+    let (a, b) = (geom_loops_store1[0][0].len(), geom_loops_store1[0][1].len());
+    let compatible0 = match geom_loops_store1[0][0].status() {
+        ShapesOpStatus::And => a == 3,
+        ShapesOpStatus::Or => a == 5,
+        _ => false,
+    };
+    let compatible1 = match geom_loops_store1[0][1].status() {
+        ShapesOpStatus::And => b == 3,
+        ShapesOpStatus::Or => b == 5,
+        _ => false,
+    };
+    assert!(
+        compatible0 && compatible1 && a * b == 15,
+        "{:?}",
+        geom_loops_store1.display(wire_id_format)
+    );
+    assert_eq!(geom_loops_store1[1].len(), 2);
+    assert!(
+        geom_loops_store1[1][0].is_closed(),
+        "{:?}",
+        geom_loops_store1.display(wire_id_format)
+    );
+    assert!(
+        geom_loops_store1[1][1].is_closed(),
+        "{:?}",
+        geom_loops_store1.display(wire_id_format)
+    );
+    assert!(
+        geom_loops_store1[1][0].is_geometric_consistent(),
+        "{:?}",
+        geom_loops_store1.display(wire_geom_format)
+    );
+    assert!(
+        geom_loops_store1[1][1].is_geometric_consistent(),
+        "{:?}",
+        geom_loops_store1.display(wire_geom_format)
+    );
+    let (a, b) = (geom_loops_store1[1][0].len(), geom_loops_store1[1][1].len());
+    let compatible0 = match geom_loops_store1[1][0].status() {
+        ShapesOpStatus::And => a == 3,
+        ShapesOpStatus::Or => a == 5,
+        _ => false,
+    };
+    let compatible1 = match geom_loops_store1[1][1].status() {
+        ShapesOpStatus::And => b == 3,
+        ShapesOpStatus::Or => b == 5,
+        _ => false,
+    };
+    assert!(
+        compatible0 && compatible1 && a * b == 15,
+        "{:?}",
+        geom_loops_store1.display(wire_id_format)
+    );
+}
+
+#[test]
+fn crossing_edges() {
+    let arc00: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new(
+        KnotVec::bezier_knot(2),
+        vec![
+            Vector4::new(1.0, 0.0, 1.0, 1.0),
+            Vector4::new(0.0, 1.0, 0.0, 0.0),
+            Vector4::new(-1.0, 0.0, 1.0, 1.0),
+        ],
+    ))
+    .into();
+    let arc01: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new(
+        KnotVec::bezier_knot(2),
+        vec![
+            Vector4::new(-1.0, 0.0, 1.0, 1.0),
+            Vector4::new(0.0, -1.0, 0.0, 0.0),
+            Vector4::new(1.0, 0.0, 1.0, 1.0),
+        ],
+    ))
+    .into();
+    let arc02: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new(
+        KnotVec::bezier_knot(2),
+        vec![
+            Vector4::new(1.0, 0.0, 1.0, 1.0),
+            Vector4::new(0.0, 0.0, -3.0, 1.0),
+            Vector4::new(-1.0, 0.0, 1.0, 1.0),
+        ],
+    ))
+    .into();
+    let arc10: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new(
+        KnotVec::bezier_knot(2),
+        vec![
+            Vector4::new(1.0, 0.0, -1.0, 1.0),
+            Vector4::new(0.0, 1.0, 0.0, 0.0),
+            Vector4::new(-1.0, 0.0, -1.0, 1.0),
+        ],
+    ))
+    .into();
+    let arc11: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new(
+        KnotVec::bezier_knot(2),
+        vec![
+            Vector4::new(-1.0, 0.0, -1.0, 1.0),
+            Vector4::new(0.0, -1.0, 0.0, 0.0),
+            Vector4::new(1.0, 0.0, -1.0, 1.0),
+        ],
+    ))
+    .into();
+    let arc12: AlternativeIntersection = NurbsCurve::new(BSplineCurve::new(
+        KnotVec::bezier_knot(2),
+        vec![
+            Vector4::new(1.0, 0.0, -1.0, 1.0),
+            Vector4::new(0.0, 0.0, 3.0, 1.0),
+            Vector4::new(-1.0, 0.0, -1.0, 1.0),
+        ],
+    ))
+    .into();
+    let (surface0, surface1) = parabola_surfaces();
+
+    let v00 = Vertex::new(Point3::new(1.0, 0.0, 1.0));
+    let v01 = Vertex::new(Point3::new(-1.0, 0.0, 1.0));
+    let edge00 = Edge::new(&v00, &v01, arc00);
+    let edge01 = Edge::new(&v01, &v00, arc01);
+    let edge02 = Edge::new(&v00, &v01, arc02);
+    let wire00: Wire<_, _> = vec![edge00, edge02.inverse()].into();
+    let wire01: Wire<_, _> = vec![edge01, edge02].into();
+    let face00 = Face::new(vec![wire00], surface0.clone());
+    let face01 = Face::new(vec![wire01], surface0);
+    let geom_shell0: Shell<_, _, _> = vec![face00.inverse(), face01.inverse()].into();
+
+    let v10 = Vertex::new(Point3::new(1.0, 0.0, -1.0));
+    let v11 = Vertex::new(Point3::new(-1.0, 0.0, -1.0));
+    let edge10 = Edge::new(&v10, &v11, arc10);
+    let edge11 = Edge::new(&v11, &v10, arc11);
+    let edge12 = Edge::new(&v10, &v11, arc12);
+    let wire10: Wire<_, _> = vec![edge10, edge12.inverse()].into();
+    let wire11: Wire<_, _> = vec![edge11, edge12].into();
+    let face10 = Face::new(vec![wire10], surface1.clone());
+    let face11 = Face::new(vec![wire11], surface1);
+    let geom_shell1: Shell<_, _, _> = vec![face10, face11].into();
+
+    let poly_shell0 = geom_shell0.triangulation(TOL);
+    let poly_shell1 = geom_shell1.triangulation(TOL);
+    let LoopsStoreQuadruple {
+        geom_loops_store0,
+        geom_loops_store1,
+        ..
+    } = create_loops_stores(&geom_shell0, &poly_shell0, &geom_shell1, &poly_shell1, TOL).unwrap();
+
+    let vertex_format = VertexDisplayFormat::AsPoint;
+    let edge_id_format = EdgeDisplayFormat::VerticesTupleAndID { vertex_format };
+    let wire_id_format = WireDisplayFormat::EdgesListTuple {
+        edge_format: edge_id_format,
+    };
+    let edge_geom_format = EdgeDisplayFormat::VerticesTupleAndCurve { vertex_format };
+    let wire_geom_format = WireDisplayFormat::EdgesListTuple {
+        edge_format: edge_geom_format,
+    };
+    assert_eq!(geom_loops_store0.len(), 2);
+    assert_eq!(geom_loops_store0[0].len(), 2);
+    assert!(
+        geom_loops_store0[0][0].is_closed(),
+        "{:?}",
+        geom_loops_store0.display(wire_id_format)
+    );
+    assert!(
+        geom_loops_store0[0][1].is_closed(),
+        "{:?}",
+        geom_loops_store0.display(wire_id_format)
+    );
+    assert!(
+        geom_loops_store0[0][0].is_geometric_consistent(),
+        "{:?}",
+        geom_loops_store0.display(wire_geom_format)
+    );
+    assert!(
+        geom_loops_store0[0][1].is_geometric_consistent(),
+        "{:?}",
+        geom_loops_store0.display(wire_geom_format)
+    );
+    let (a, b) = (geom_loops_store0[0][0].len(), geom_loops_store0[0][1].len());
+    let compatible0 = match geom_loops_store0[0][0].status() {
+        ShapesOpStatus::And => a == 2,
+        ShapesOpStatus::Or => a == 4,
+        _ => false,
+    };
+    let compatible1 = match geom_loops_store0[0][1].status() {
+        ShapesOpStatus::And => b == 2,
+        ShapesOpStatus::Or => b == 4,
+        _ => false,
+    };
+    assert!(
+        compatible0 && compatible1 && a * b == 8,
+        "{:?}",
+        geom_loops_store0.display(wire_id_format)
+    );
+    assert_eq!(geom_loops_store0[1].len(), 2);
+    assert!(
+        geom_loops_store0[1][0].is_closed(),
+        "{:?}",
+        geom_loops_store0.display(wire_id_format)
+    );
+    assert!(
+        geom_loops_store0[1][1].is_closed(),
+        "{:?}",
+        geom_loops_store0.display(wire_id_format)
+    );
+    assert!(
+        geom_loops_store0[1][0].is_geometric_consistent(),
+        "{:?}",
+        geom_loops_store0.display(wire_geom_format)
+    );
+    assert!(
+        geom_loops_store0[1][1].is_geometric_consistent(),
+        "{:?}",
+        geom_loops_store0.display(wire_geom_format)
+    );
+    let (a, b) = (geom_loops_store0[1][0].len(), geom_loops_store0[1][1].len());
+    let compatible0 = match geom_loops_store0[1][0].status() {
+        ShapesOpStatus::And => a == 2,
+        ShapesOpStatus::Or => a == 4,
+        _ => false,
+    };
+    let compatible1 = match geom_loops_store0[1][1].status() {
+        ShapesOpStatus::And => b == 2,
+        ShapesOpStatus::Or => b == 4,
+        _ => false,
+    };
+    assert!(
+        compatible0 && compatible1 && a * b == 8,
+        "{:?}",
+        geom_loops_store0.display(wire_id_format)
+    );
+    assert_eq!(geom_loops_store1.len(), 2);
+    assert_eq!(geom_loops_store1[0].len(), 2);
+    assert!(
+        geom_loops_store1[0][0].is_closed(),
+        "{:?}",
+        geom_loops_store1.display(wire_id_format)
+    );
+    assert!(
+        geom_loops_store1[0][1].is_closed(),
+        "{:?}",
+        geom_loops_store1.display(wire_id_format)
+    );
+    assert!(
+        geom_loops_store1[0][0].is_geometric_consistent(),
+        "{:?}",
+        geom_loops_store1.display(wire_geom_format)
+    );
+    assert!(
+        geom_loops_store1[0][1].is_geometric_consistent(),
+        "{:?}",
+        geom_loops_store1.display(wire_geom_format)
+    );
+    let (a, b) = (geom_loops_store1[0][0].len(), geom_loops_store1[0][1].len());
+    let compatible0 = match geom_loops_store1[0][0].status() {
+        ShapesOpStatus::And => a == 2,
+        ShapesOpStatus::Or => a == 4,
+        _ => false,
+    };
+    let compatible1 = match geom_loops_store1[0][1].status() {
+        ShapesOpStatus::And => b == 2,
+        ShapesOpStatus::Or => b == 4,
+        _ => false,
+    };
+    assert!(
+        compatible0 && compatible1 && a * b == 8,
+        "{:?}",
+        geom_loops_store1.display(wire_id_format)
+    );
+    assert_eq!(geom_loops_store1[1].len(), 2);
+    assert!(
+        geom_loops_store1[1][0].is_closed(),
+        "{:?}",
+        geom_loops_store1.display(wire_id_format)
+    );
+    assert!(
+        geom_loops_store1[1][1].is_closed(),
+        "{:?}",
+        geom_loops_store1.display(wire_id_format)
+    );
+    assert!(
+        geom_loops_store1[1][0].is_geometric_consistent(),
+        "{:?}",
+        geom_loops_store1.display(wire_geom_format)
+    );
+    assert!(
+        geom_loops_store1[1][1].is_geometric_consistent(),
+        "{:?}",
+        geom_loops_store1.display(wire_geom_format)
+    );
+    let (a, b) = (geom_loops_store1[1][0].len(), geom_loops_store1[1][1].len());
+    let compatible0 = match geom_loops_store1[1][0].status() {
+        ShapesOpStatus::And => a == 2,
+        ShapesOpStatus::Or => a == 4,
+        _ => false,
+    };
+    let compatible1 = match geom_loops_store1[1][1].status() {
+        ShapesOpStatus::And => b == 2,
+        ShapesOpStatus::Or => b == 4,
+        _ => false,
+    };
+    assert!(
+        compatible0 && compatible1 && a * b == 8,
+        "{:?}",
+        geom_loops_store0.display(wire_id_format)
+    );
+}
@@ -0,0 +1,7 @@
+mod divide_face;
+mod faces_classification;
+mod integrate;
+mod intersection_curve;
+mod loops_store;
+mod polyline_construction;
+pub use integrate::{and, or, ShapeOpsCurve, ShapeOpsSurface};
@@ -0,0 +1,117 @@
+use rustc_hash::{FxHashMap as HashMap, FxHashSet as HashSet};
+use std::collections::VecDeque;
+use truck_base::{cgmath64::*, tolerance::*};
+use truck_meshalgo::prelude::PolylineCurve;
+
+pub fn construct_polylines(lines: &[(Point3, Point3)]) -> Vec<PolylineCurve<Point3>> {
+    let mut graph: Graph = lines.iter().collect();
+    let mut res = Vec::new();
+    while !graph.is_empty() {
+        let (mut idx, node) = graph.get_one();
+        let mut wire: VecDeque<_> = vec![node.coord].into();
+        while let Some((idx0, pt)) = graph.get_a_next_node(idx) {
+            idx = idx0;
+            wire.push_back(pt);
+        }
+        let mut idx = PointIndex::from(wire[0]);
+        while let Some((idx0, pt)) = graph.get_a_next_node(idx) {
+            idx = idx0;
+            wire.push_front(pt);
+        }
+        res.push(PolylineCurve(wire.into()));
+    }
+    res
+}
+
+#[derive(Clone, Debug, Copy, Hash, PartialEq, Eq)]
+struct PointIndex([i64; 3]);
+
+impl From<Point3> for PointIndex {
+    #[inline(always)]
+    fn from(pt: Point3) -> PointIndex {
+        let idx = pt.add_element_wise(TOLERANCE) / (2.0 * TOLERANCE);
+        PointIndex(idx.cast::<i64>().unwrap().into())
+    }
+}
+
+struct Node {
+    coord: Point3,
+    adjacency: HashSet<PointIndex>,
+}
+
+impl Node {
+    #[inline(always)]
+    fn new(coord: Point3, adjacency: HashSet<PointIndex>) -> Node { Node { coord, adjacency } }
+
+    fn pop_one_adjacency(&mut self) -> PointIndex {
+        let idx = *self.adjacency.iter().next().unwrap();
+        self.adjacency.remove(&idx);
+        idx
+    }
+}
+
+struct Graph(HashMap<PointIndex, Node>);
+
+impl std::ops::Deref for Graph {
+    type Target = HashMap<PointIndex, Node>;
+    #[inline(always)]
+    fn deref(&self) -> &Self::Target { &self.0 }
+}
+
+impl std::ops::DerefMut for Graph {
+    #[inline(always)]
+    fn deref_mut(&mut self) -> &mut Self::Target { &mut self.0 }
+}
+
+impl Graph {
+    fn add_half_edge(&mut self, pt0: Point3, pt1: Point3) {
+        let idx0 = pt0.into();
+        let idx1 = pt1.into();
+        if let Some(node) = self.get_mut(&idx0) {
+            node.adjacency.insert(idx1);
+        } else {
+            let mut set = HashSet::default();
+            set.insert(idx1);
+            self.insert(idx0, Node::new(pt0, set));
+        }
+    }
+
+    fn add_edge(&mut self, line: (Point3, Point3)) {
+        if !line.0.near(&line.1) {
+            self.add_half_edge(line.0, line.1);
+            self.add_half_edge(line.1, line.0);
+        }
+    }
+
+    #[inline(always)]
+    fn get_one(&self) -> (PointIndex, &Node) {
+        let (idx, node) = self.iter().next().unwrap();
+        (*idx, node)
+    }
+
+    fn get_a_next_node(&mut self, idx: PointIndex) -> Option<(PointIndex, Point3)> {
+        let node = self.get_mut(&idx)?;
+        let idx0 = node.pop_one_adjacency();
+        if node.adjacency.is_empty() {
+            self.remove(&idx);
+        }
+        let node = self.get_mut(&idx0)?;
+        node.adjacency.remove(&idx);
+        let pt = node.coord;
+        if node.adjacency.is_empty() {
+            self.remove(&idx0);
+        }
+        Some((idx0, pt))
+    }
+}
+
+impl<'a> FromIterator<&'a (Point3, Point3)> for Graph {
+    fn from_iter<I: IntoIterator<Item = &'a (Point3, Point3)>>(iter: I) -> Graph {
+        let mut res = Graph(HashMap::default());
+        iter.into_iter().for_each(|line| res.add_edge(*line));
+        res
+    }
+}
+
+#[cfg(test)]
+mod tests;
@@ -0,0 +1,113 @@
+use super::*;
+
+#[test]
+fn construct_polylines_positive0() {
+    let lines = vec![
+        (Point3::new(0.0, 0.0, 0.0), Point3::new(1.0, 0.0, 0.0)),
+        (Point3::new(1.0, 0.0, 0.0), Point3::new(0.0, 1.0, 0.0)),
+        (Point3::new(1.0, 1.0, 0.0), Point3::new(0.0, 0.0, 1.0)),
+        (Point3::new(0.0, 1.0, 1.0), Point3::new(1.0, 1.0, 1.0)),
+        (Point3::new(0.0, 0.0, 1.0), Point3::new(1.0, 0.0, 1.0)),
+        (Point3::new(0.0, 1.0, 0.0), Point3::new(1.0, 1.0, 0.0)),
+        (Point3::new(1.0, 1.0, 1.0), Point3::new(0.0, 0.0, 0.0)),
+        (Point3::new(1.0, 0.0, 1.0), Point3::new(0.0, 1.0, 1.0)),
+    ];
+    let polyline = construct_polylines(&lines);
+    assert_eq!(polyline.len(), 1);
+    assert_eq!(polyline[0].len(), 9);
+
+    let mut sign = None;
+    for line in polyline[0].windows(2) {
+        let a = line[0][0] + line[0][1] * 2.0 + line[0][2] * 4.0;
+        let b = line[1][0] + line[1][1] * 2.0 + line[1][2] * 4.0;
+        let x = b - a;
+        assert!(f64::abs(x) == 1.0 || f64::abs(x) == 7.0);
+        let s = f64::signum(x * (x - 2.0) * (x + 2.0));
+        if let Some(sign) = sign {
+            assert!(s == sign);
+        } else {
+            sign = Some(s);
+        }
+    }
+}
+
+#[test]
+fn construct_polylines_positive1() {
+    let lines = vec![
+        (Point3::new(0.0, 0.0, 0.0), Point3::new(1.0, 0.0, 0.0)),
+        (Point3::new(1.0, 0.0, 0.0), Point3::new(1.0, 1.0, 0.0)),
+        (Point3::new(1.0, 0.0, 1.0), Point3::new(1.0, 1.0, 1.0)),
+        (Point3::new(1.0, 1.0, 1.0), Point3::new(0.0, 1.0, 1.0)),
+        (Point3::new(1.0, 1.0, 0.0), Point3::new(0.0, 1.0, 0.0)),
+        (Point3::new(0.0, 1.0, 0.0), Point3::new(0.0, 0.0, 0.0)),
+        (Point3::new(0.0, 0.0, 1.0), Point3::new(1.0, 0.0, 1.0)),
+        (Point3::new(0.0, 1.0, 1.0), Point3::new(0.0, 0.0, 1.0)),
+    ];
+    let polyline = construct_polylines(&lines);
+    assert_eq!(polyline.len(), 2);
+    assert_eq!(polyline[0].len(), 5);
+    assert_eq!(polyline[1].len(), 5);
+}
+
+#[test]
+fn construct_polylines_positive2() {
+    let lines = vec![
+        (Point3::new(0.0, 0.0, 0.0), Point3::new(1.0, 0.0, 0.0)),
+        (Point3::new(1.0, 0.0, 0.0), Point3::new(0.0, 1.0, 0.0)),
+        (Point3::new(1.0, 1.0, 0.0), Point3::new(0.0, 0.0, 1.0)),
+        (Point3::new(0.0, 1.0, 1.0), Point3::new(1.0, 1.0, 1.0)),
+        (Point3::new(0.0, 0.0, 1.0), Point3::new(1.0, 0.0, 1.0)),
+        (Point3::new(1.0, 1.0, 0.0), Point3::new(1.0, 1.0, 0.0)),
+        (Point3::new(0.0, 1.0, 0.0), Point3::new(1.0, 1.0, 0.0)),
+        (Point3::new(1.0, 1.0, 1.0), Point3::new(0.0, 0.0, 0.0)),
+        (Point3::new(1.0, 0.0, 1.0), Point3::new(0.0, 1.0, 1.0)),
+    ];
+    let polyline = construct_polylines(&lines);
+    assert_eq!(polyline.len(), 1);
+    assert_eq!(polyline[0].len(), 9);
+
+    let mut sign = None;
+    for line in polyline[0].windows(2) {
+        let a = line[0][0] + line[0][1] * 2.0 + line[0][2] * 4.0;
+        let b = line[1][0] + line[1][1] * 2.0 + line[1][2] * 4.0;
+        let x = b - a;
+        assert!(f64::abs(x) == 1.0 || f64::abs(x) == 7.0);
+        let s = f64::signum(x * (x - 2.0) * (x + 2.0));
+        if let Some(sign) = sign {
+            assert!(s == sign);
+        } else {
+            sign = Some(s);
+        }
+    }
+}
+
+#[test]
+fn construct_polylines_positive3() {
+    let lines = vec![
+        (Point3::new(0.0, 0.0, 0.0), Point3::new(1.0, 0.0, 0.0)),
+        (Point3::new(1.0, 0.0, 0.0), Point3::new(0.0, 1.0, 0.0)),
+        (Point3::new(1.0, 1.0, 0.0), Point3::new(0.0, 0.0, 1.0)),
+        (Point3::new(0.0, 1.0, 1.0), Point3::new(1.0, 1.0, 1.0)),
+        (Point3::new(0.0, 0.0, 1.0), Point3::new(1.0, 0.0, 1.0)),
+        (Point3::new(1.0, 1.0, 0.0), Point3::new(1.0, 1.0, 0.0)),
+        (Point3::new(0.0, 1.0, 0.0), Point3::new(1.0, 1.0, 0.0)),
+        (Point3::new(1.0, 0.0, 1.0), Point3::new(0.0, 1.0, 1.0)),
+    ];
+    let polyline = construct_polylines(&lines);
+    assert_eq!(polyline.len(), 1);
+    assert_eq!(polyline[0].len(), 8);
+
+    let mut sign = None;
+    for line in polyline[0].windows(2) {
+        let a = line[0][0] + line[0][1] * 2.0 + line[0][2] * 4.0;
+        let b = line[1][0] + line[1][1] * 2.0 + line[1][2] * 4.0;
+        let x = b - a;
+        assert!(f64::abs(x) == 1.0);
+        let s = f64::signum(x * (x - 2.0) * (x + 2.0));
+        if let Some(sign) = sign {
+            assert!(s == sign);
+        } else {
+            sign = Some(s);
+        }
+    }
+}