Merge pull request #1209 from bernhard-42/intersect-optimized

Intersect optimized

src/build123d/geometry.py

@@ -1182,6 +1182,20 @@ class BoundBox:
             and second_box.max.Z < self.max.Z
         )
 
+    def overlaps(self, other: BoundBox, tolerance: float = TOLERANCE) -> bool:
+        """Check if this bounding box overlaps with another.
+
+        Args:
+            other: BoundBox to check overlap with
+            tolerance: Distance tolerance for overlap detection
+
+        Returns:
+            True if bounding boxes overlap (share any volume), False otherwise
+        """
+        if self.wrapped is None or other.wrapped is None:
+            return False
+        return self.wrapped.Distance(other.wrapped) <= tolerance
+
     def to_align_offset(self, align: Align2DType | Align3DType) -> Vector:
         """Amount to move object to achieve the desired alignment"""
         return to_align_offset(self.min, self.max, align)

src/build123d/operations_part.py

@@ -230,16 +230,14 @@ def extrude(
                     )
                 )
 
+    if both and len(new_solids) > 1:
+        fused_solids = new_solids.pop().fuse(*new_solids)
+        new_solids = fused_solids if isinstance(fused_solids, list) else [fused_solids]
+    if clean:
+        new_solids = [solid.clean() for solid in new_solids]
+
     if context is not None:
         context._add_to_context(*new_solids, clean=clean, mode=mode)
-    else:
-        if len(new_solids) > 1:
-            fused_solids = new_solids.pop().fuse(*new_solids)
-            new_solids = (
-                fused_solids if isinstance(fused_solids, list) else [fused_solids]
-            )
-        if clean:
-            new_solids = [solid.clean() for solid in new_solids]
 
     return Part(ShapeList(new_solids).solids())
 

src/build123d/topology/composite.py

@@ -713,149 +713,88 @@ class Compound(Mixin3D[TopoDS_Compound]):
 
         return results
 
-    def intersect(
-        self, *to_intersect: Shape | Vector | Location | Axis | Plane
-    ) -> None | ShapeList[Vertex | Edge | Face | Solid]:
-        """Intersect Compound with Shape or geometry object
+    def _intersect(
+        self,
+        other: Shape | Vector | Location | Axis | Plane,
+        tolerance: float = 1e-6,
+        include_touched: bool = False,
+    ) -> ShapeList | None:
+        """Single-object intersection for Compound (OR semantics).
+
+        Distributes intersection over elements, collecting all results:
+            Compound([a, b]).intersect(s) = (a ∩ s) ∪ (b ∩ s)
+            Compound([a, b]).intersect(Compound([c, d])) = (a ∩ c) ∪ (a ∩ d) ∪ (b ∩ c) ∪ (b ∩ d)
+
+        Handles both build123d assemblies (children) and OCCT Compounds (list()).
+        Nested Compounds are handled by recursion.
 
         Args:
-            to_intersect (Shape | Vector | Location | Axis | Plane): objects to intersect
+            other: Shape or geometry object to intersect with
+            tolerance: tolerance for intersection detection
+            include_touched: if True, include boundary contacts
+                (only relevant when Solids are involved)
+        """
+        # Convert geometry objects
+        if isinstance(other, Vector):
+            other = Vertex(other)
+        elif isinstance(other, Location):
+            other = Vertex(other.position)
+        elif isinstance(other, Axis):
+            other = Edge(other)
+        elif isinstance(other, Plane):
+            other = Face(other)
+
+        # Get self elements: assembly children or OCCT direct children
+        self_elements = self.children if self.children else list(self)
+
+        if not self_elements:
+            return None
+
+        results: ShapeList = ShapeList()
+
+        # Distribute over elements (OR semantics for Compound arguments)
+        if isinstance(other, Compound):
+            other_elements = other.children if other.children else list(other)
+        else:
+            other_elements = [other]
+
+        for self_elem in self_elements:
+            for other_elem in other_elements:
+                intersection = self_elem._intersect(
+                    other_elem, tolerance, include_touched
+                )
+                if intersection:
+                    results.extend(intersection)
+
+        # Remove duplicates using Shape's __hash__
+        unique = ShapeList(set(results))
+
+        return unique if unique else None
+
+    def touch(
+        self, other: Shape, tolerance: float = 1e-6
+    ) -> ShapeList[Vertex | Edge | Face]:
+        """Distribute touch over compound elements.
+
+        Iterates over elements and collects touch results. Only Solid and
+        Face elements produce boundary contacts; other shapes return empty.
+
+        Args:
+            other: Shape to check boundary contacts with
+            tolerance: tolerance for contact detection
 
         Returns:
-            ShapeList[Vertex | Edge | Face | Solid] | None: ShapeList of vertices, edges,
-                faces, and/or solids.
+            ShapeList of boundary contact geometry (empty if no contact)
         """
+        results: ShapeList = ShapeList()
 
-        def to_vector(objs: Iterable) -> ShapeList:
-            return ShapeList([Vector(v) if isinstance(v, Vertex) else v for v in objs])
+        # Get elements: assembly children or OCCT direct children
+        elements = self.children if self.children else list(self)
 
-        def to_vertex(objs: Iterable) -> ShapeList:
-            return ShapeList([Vertex(v) if isinstance(v, Vector) else v for v in objs])
+        for elem in elements:
+            results.extend(elem.touch(other, tolerance))
 
-        def bool_op(
-            args: Sequence,
-            tools: Sequence,
-            operation: BRepAlgoAPI_Common | BRepAlgoAPI_Section,
-        ) -> ShapeList:
-            # Wrap Shape._bool_op for corrected output
-            intersections: Shape | ShapeList = Shape()._bool_op(args, tools, operation)
-            if isinstance(intersections, ShapeList):
-                return intersections
-            if isinstance(intersections, Shape) and not intersections.is_null:
-                return ShapeList([intersections])
-            return ShapeList()
-
-        def expand_compound(compound: Compound) -> ShapeList:
-            shapes = ShapeList(compound.children)
-            for shape_type in [Vertex, Edge, Wire, Face, Shell, Solid]:
-                shapes.extend(compound.get_type(shape_type))
-            return shapes
-
-        def filter_shapes_by_order(shapes: ShapeList, orders: list) -> ShapeList:
-            # Remove lower order shapes from list which *appear* to be part of
-            # a higher order shape using a lazy distance check
-            # (sufficient for vertices, may be an issue for higher orders)
-            order_groups = []
-            for order in orders:
-                order_groups.append(
-                    ShapeList([s for s in shapes if isinstance(s, order)])
-                )
-
-            filtered_shapes = order_groups[-1]
-            for i in range(len(order_groups) - 1):
-                los = order_groups[i]
-                his: list = sum(order_groups[i + 1 :], [])
-                filtered_shapes.extend(
-                    ShapeList(
-                        lo
-                        for lo in los
-                        if all(lo.distance_to(hi) > TOLERANCE for hi in his)
-                    )
-                )
-
-            return filtered_shapes
-
-        common_set: ShapeList[Shape] = expand_compound(self)
-        target: ShapeList | Shape
-        for other in to_intersect:
-            # Conform target type
-            match other:
-                case Axis():
-                    # BRepAlgoAPI_Section seems happier if Edge isnt infinite
-                    bbox = self.bounding_box()
-                    dist = self.distance_to(other.position)
-                    dist = dist if dist >= 1 else 1
-                    target = Edge.make_line(
-                        other.position - other.direction * bbox.diagonal * dist,
-                        other.position + other.direction * bbox.diagonal * dist,
-                    )
-                case Plane():
-                    target = Face(other)
-                case Vector():
-                    target = Vertex(other)
-                case Location():
-                    target = Vertex(other.position)
-                case Compound():
-                    target = expand_compound(other)
-                case _ if issubclass(type(other), Shape):
-                    target = other
-                case _:
-                    raise ValueError(f"Unsupported type to_intersect: {type(other)}")
-
-            # Find common matches
-            common: list[Vertex | Edge | Wire | Face | Shell | Solid] = []
-            result: ShapeList
-            for obj in common_set:
-                if isinstance(target, Shape):
-                    target = ShapeList([target])
-                result = ShapeList()
-                for t in target:
-                    operation: BRepAlgoAPI_Section | BRepAlgoAPI_Common = (
-                        BRepAlgoAPI_Section()
-                    )
-                    result.extend(bool_op((obj,), (t,), operation))
-                    if (
-                        not isinstance(obj, Edge | Wire)
-                        and not isinstance(t, Edge | Wire)
-                    ) or (
-                        isinstance(obj, Solid | Compound)
-                        or isinstance(t, Solid | Compound)
-                    ):
-                        # Face + Edge combinations may produce an intersection
-                        # with Common but always with Section.
-                        # No easy way to deduplicate
-                        # Many Solid + Edge combinations need Common
-                        operation = BRepAlgoAPI_Common()
-                        result.extend(bool_op((obj,), (t,), operation))
-
-                if result:
-                    common.extend(result)
-
-            expanded: ShapeList = ShapeList()
-            if common:
-                for shape in common:
-                    if isinstance(shape, Compound):
-                        expanded.extend(expand_compound(shape))
-                    else:
-                        expanded.append(shape)
-
-            if expanded:
-                common_set = ShapeList()
-                for shape in expanded:
-                    if isinstance(shape, Wire):
-                        common_set.extend(shape.edges())
-                    elif isinstance(shape, Shell):
-                        common_set.extend(shape.faces())
-                    else:
-                        common_set.append(shape)
-                common_set = to_vertex(set(to_vector(common_set)))
-                common_set = filter_shapes_by_order(
-                    common_set, [Vertex, Edge, Face, Solid]
-                )
-            else:
-                return None
-
-        return ShapeList(common_set)
+        return ShapeList(set(results))
 
     def project_to_viewport(
         self,

src/build123d/topology/one_d.py

@@ -713,126 +713,73 @@ class Mixin1D(Shape[TOPODS]):
 
         return Vector(curve.Value(umax))
 
-    def intersect(
-        self, *to_intersect: Shape | Vector | Location | Axis | Plane
-    ) -> None | ShapeList[Vertex | Edge]:
-        """Intersect Edge with Shape or geometry object
+    def _intersect(
+        self,
+        other: Shape | Vector | Location | Axis | Plane,
+        tolerance: float = 1e-6,
+        include_touched: bool = False,
+    ) -> ShapeList | None:
+        """Single-object intersection for Edge/Wire.
+
+        Returns same-dimension overlap or crossing geometry:
+        - 1D + 1D → Edge (collinear overlap) + Vertex (crossing)
+        - 1D + Face/Solid/Compound → delegates to other._intersect(self)
 
         Args:
-            to_intersect (Shape | Vector | Location | Axis | Plane): objects to intersect
-
-        Returns:
-            ShapeList[Vertex | Edge] | None: ShapeList of vertices and/or edges
+            other: Shape or geometry object to intersect with
+            tolerance: tolerance for intersection detection
+            include_touched: if True, include boundary contacts
+                (only relevant when Solids are involved)
         """
+        # Convert geometry objects to shapes
+        if isinstance(other, Vector):
+            other = Vertex(other)
+        elif isinstance(other, Location):
+            other = Vertex(other.position)
+        elif isinstance(other, Axis):
+            other = Edge(other)
 
-        def to_vector(objs: Iterable) -> ShapeList:
-            return ShapeList([Vector(v) if isinstance(v, Vertex) else v for v in objs])
+        results: ShapeList = ShapeList()
 
-        def to_vertex(objs: Iterable) -> ShapeList:
-            return ShapeList([Vertex(v) if isinstance(v, Vector) else v for v in objs])
+        # Trim infinite edges before OCCT operations
+        if isinstance(other, Edge) and other.is_infinite:
+            bbox = self.bounding_box(optimal=False)
+            other = other.trim_infinite(
+                bbox.diagonal + (other.center() - bbox.center()).length
+            )
 
-        def bool_op(
-            args: Sequence,
-            tools: Sequence,
-            operation: BRepAlgoAPI_Section | BRepAlgoAPI_Common,
-        ) -> ShapeList:
-            # Wrap Shape._bool_op for corrected output
-            intersections: Shape | ShapeList = Shape()._bool_op(args, tools, operation)
-            if isinstance(intersections, ShapeList):
-                return intersections or ShapeList()
-            if isinstance(intersections, Shape) and not intersections.is_null:
-                return ShapeList([intersections])
-            return ShapeList()
+        # 1D + Plane: run Section directly with OCP Face
+        if isinstance(other, Plane):
+            face: Shape = Shape(BRepBuilderAPI_MakeFace(other.wrapped).Face())
+            section = self._bool_op_list((self,), (face,), BRepAlgoAPI_Section())
+            results.extend(section.expand())
 
-        def filter_shapes_by_order(shapes: ShapeList, orders: list) -> ShapeList:
-            # Remove lower order shapes from list which *appear* to be part of
-            # a higher order shape using a lazy distance check
-            # (sufficient for vertices, may be an issue for higher orders)
-            order_groups = []
-            for order in orders:
-                order_groups.append(
-                    ShapeList([s for s in shapes if isinstance(s, order)])
-                )
+        # 1D + 1D: Common (collinear overlap) + Section (crossing vertices)
+        elif isinstance(other, (Edge, Wire)):
+            common = self._bool_op_list(
+                (self,), (other,), BRepAlgoAPI_Common()
+            )
+            results.extend(common.expand())
+            section = self._bool_op_list(
+                (self,), (other,), BRepAlgoAPI_Section()
+            )
+            # Extract vertices from section (edges already in Common for wires)
+            for shape in section:
+                if isinstance(shape, Vertex) and not shape.is_null:
+                    results.append(shape)
 
-            filtered_shapes = order_groups[-1]
-            for i in range(len(order_groups) - 1):
-                los = order_groups[i]
-                his: list = sum(order_groups[i + 1 :], [])
-                filtered_shapes.extend(
-                    ShapeList(
-                        lo
-                        for lo in los
-                        if all(lo.distance_to(hi) > TOLERANCE for hi in his)
-                    )
-                )
+        # 1D + Vertex: point containment on edge
+        elif isinstance(other, Vertex):
+            if other.distance_to(self) <= tolerance:
+                results.append(other)
 
-            return filtered_shapes
+        # Delegate to higher-order shapes (Face, Solid, etc.)
+        else:
+            result = other._intersect(self, tolerance, include_touched)
+            if result:
+                results.extend(result)
 
-        common_set: ShapeList[Vertex | Edge | Wire] = ShapeList([self])
-        target: Shape | Plane
-        for other in to_intersect:
-            # Conform target type
-            match other:
-                case Axis():
-                    # BRepAlgoAPI_Section seems happier if Edge isnt infinite
-                    bbox = self.bounding_box()
-                    dist = self.distance_to(other.position)
-                    dist = dist if dist >= 1 else 1
-                    target = Edge.make_line(
-                        other.position - other.direction * bbox.diagonal * dist,
-                        other.position + other.direction * bbox.diagonal * dist,
-                    )
-                case Plane():
-                    target = other
-                case Vector():
-                    target = Vertex(other)
-                case Location():
-                    target = Vertex(other.position)
-                case _ if issubclass(type(other), Shape):
-                    target = other
-                case _:
-                    raise ValueError(f"Unsupported type to_intersect: {type(other)}")
-
-            # Find common matches
-            common: list[Vertex | Edge | Wire] = []
-            result: ShapeList | None
-            for obj in common_set:
-                match (obj, target):
-                    case (_, Plane()):
-                        assert isinstance(other.wrapped, gp_Pln)
-                        target = Shape(BRepBuilderAPI_MakeFace(other.wrapped).Face())
-                        operation1 = BRepAlgoAPI_Section()
-                        result = bool_op((obj,), (target,), operation1)
-                        operation2 = BRepAlgoAPI_Common()
-                        result.extend(bool_op((obj,), (target,), operation2))
-
-                    case (_, Vertex() | Edge() | Wire()):
-                        operation1 = BRepAlgoAPI_Section()
-                        section = bool_op((obj,), (target,), operation1)
-                        result = section
-                        if not section:
-                            operation2 = BRepAlgoAPI_Common()
-                            result.extend(bool_op((obj,), (target,), operation2))
-
-                    case _ if issubclass(type(target), Shape):
-                        result = target.intersect(obj)
-
-                if result:
-                    common.extend(result)
-
-            if common:
-                common_set = ShapeList()
-                for shape in common:
-                    if isinstance(shape, Wire):
-                        common_set.extend(shape.edges())
-                    else:
-                        common_set.append(shape)
-                common_set = to_vertex(set(to_vector(common_set)))
-                common_set = filter_shapes_by_order(common_set, [Vertex, Edge])
-            else:
-                return None
-
-        return ShapeList(common_set)
+        return results if results else None
 
     def location_at(
         self,
@@ -2692,6 +2639,132 @@ class Edge(Mixin1D[TopoDS_Edge]):
             raise ValueError("Can't find adaptor for empty edge")
         return BRepAdaptor_Curve(self.wrapped)
 
+    def geom_equal(
+        self,
+        other: Edge,
+        tol: float = 1e-6,
+        num_interpolation_points: int = 5,
+    ) -> bool:
+        """Compare two edges for geometric equality within tolerance.
+
+        This compares the geometric properties of two edges, not their topological
+        identity. Two independently created edges with the same geometry will
+        return True.
+
+        Args:
+            other: Edge to compare with
+            tol: Tolerance for numeric comparisons. Defaults to 1e-6.
+            num_interpolation_points: Number of points to sample for unknown
+                curve types. Defaults to 5.
+
+        Returns:
+            bool: True if edges are geometrically equal within tolerance
+        """
+        if not isinstance(other, Edge):
+            return False
+
+        # geom_type must match
+        if self.geom_type != other.geom_type:
+            return False
+
+        # Common: start and end points
+        if (self @ 0) != (other @ 0) or (self @ 1) != (other @ 1):
+            return False
+
+        ga1 = self.geom_adaptor()
+        ga2 = other.geom_adaptor()
+
+        match self.geom_type:
+            case GeomType.LINE:
+                # Line: fully defined by endpoints (already checked)
+                return True
+
+            case GeomType.CIRCLE:
+                c1, c2 = ga1.Circle(), ga2.Circle()
+                return (
+                    abs(c1.Radius() - c2.Radius()) < tol
+                    and Vector(c1.Location()) == Vector(c2.Location())
+                    and Vector(c1.Axis().Direction()) == Vector(c2.Axis().Direction())
+                )
+
+            case GeomType.ELLIPSE:
+                e1, e2 = ga1.Ellipse(), ga2.Ellipse()
+                return (
+                    abs(e1.MajorRadius() - e2.MajorRadius()) < tol
+                    and abs(e1.MinorRadius() - e2.MinorRadius()) < tol
+                    and Vector(e1.Location()) == Vector(e2.Location())
+                    and Vector(e1.Axis().Direction()) == Vector(e2.Axis().Direction())
+                )
+
+            case GeomType.HYPERBOLA:
+                h1, h2 = ga1.Hyperbola(), ga2.Hyperbola()
+                return (
+                    abs(h1.MajorRadius() - h2.MajorRadius()) < tol
+                    and abs(h1.MinorRadius() - h2.MinorRadius()) < tol
+                    and Vector(h1.Location()) == Vector(h2.Location())
+                    and Vector(h1.Axis().Direction()) == Vector(h2.Axis().Direction())
+                )
+
+            case GeomType.PARABOLA:
+                p1, p2 = ga1.Parabola(), ga2.Parabola()
+                return (
+                    abs(p1.Focal() - p2.Focal()) < tol
+                    and Vector(p1.Location()) == Vector(p2.Location())
+                    and Vector(p1.Axis().Direction()) == Vector(p2.Axis().Direction())
+                )
+
+            case GeomType.BEZIER:
+                b1, b2 = ga1.Bezier(), ga2.Bezier()
+                if b1.Degree() != b2.Degree() or b1.NbPoles() != b2.NbPoles():
+                    return False
+                for i in range(1, b1.NbPoles() + 1):
+                    if Vector(b1.Pole(i)) != Vector(b2.Pole(i)):
+                        return False
+                    if b1.IsRational() and abs(b1.Weight(i) - b2.Weight(i)) >= tol:
+                        return False
+                return True
+
+            case GeomType.BSPLINE:
+                s1, s2 = ga1.BSpline(), ga2.BSpline()
+                if s1.Degree() != s2.Degree():
+                    return False
+                if s1.IsPeriodic() != s2.IsPeriodic():
+                    return False
+                if s1.NbPoles() != s2.NbPoles() or s1.NbKnots() != s2.NbKnots():
+                    return False
+                for i in range(1, s1.NbPoles() + 1):
+                    if Vector(s1.Pole(i)) != Vector(s2.Pole(i)):
+                        return False
+                    if s1.IsRational() and abs(s1.Weight(i) - s2.Weight(i)) >= tol:
+                        return False
+                for i in range(1, s1.NbKnots() + 1):
+                    if abs(s1.Knot(i) - s2.Knot(i)) >= tol:
+                        return False
+                    if s1.Multiplicity(i) != s2.Multiplicity(i):
+                        return False
+                return True
+
+            case GeomType.OFFSET:
+                oc1, oc2 = ga1.OffsetCurve(), ga2.OffsetCurve()
+                # Compare offset values and directions
+                if abs(oc1.Offset() - oc2.Offset()) >= tol:
+                    return False
+                if Vector(oc1.Direction()) != Vector(oc2.Direction()):
+                    return False
+                # Compare basis curves (recursive)
+                basis1 = Edge(BRepBuilderAPI_MakeEdge(oc1.BasisCurve()).Edge())
+                basis2 = Edge(BRepBuilderAPI_MakeEdge(oc2.BasisCurve()).Edge())
+                return basis1.geom_equal(basis2, tol)
+
+            case _:  # pragma: no cover
+                # I don't think, GeomAbs_OtherCurve can be created in Python
+                # OTHER/unknown: compare sample points
+                for i in range(1, num_interpolation_points + 1):
+                    t = i / (num_interpolation_points + 1)
+                    if (self @ t) != (other @ t):
+                        return False
+                return True
+
     def _occt_param_at(
         self, position: float, position_mode: PositionMode = PositionMode.PARAMETER
     ) -> tuple[BRepAdaptor_Curve, float, bool]:
@@ -3051,6 +3124,36 @@ class Edge(Mixin1D[TopoDS_Edge]):
         new_edge = BRepBuilderAPI_MakeEdge(trimmed_curve).Edge()
         return Edge(new_edge)
 
+    @property
+    def is_infinite(self) -> bool:
+        """Check if edge is infinite (LINE with length > 1e100)."""
+        return self.geom_type == GeomType.LINE and self.length > 1e100
+
+    def trim_infinite(self, half_length: float) -> Edge:
+        """Trim an infinite line edge to a finite length.
+
+        OCCT's boolean operations struggle with very long edges (length > 1e100).
+        This method trims such edges to a reasonable size centered at edge.center().
+
+        For non-infinite edges, returns self unchanged.
+
+        Args:
+            half_length: Half-length of the resulting edge
+
+        Returns:
+            Trimmed edge if infinite, otherwise self
+        """
+        if not self.is_infinite:
+            return self
+
+        origin = self.center()
+        direction = (self.end_point() - self.start_point()).normalized()
+
+        return Edge.make_line(
+            origin - direction * half_length,
+            origin + direction * half_length,
+        )
+
 
 class Wire(Mixin1D[TopoDS_Wire]):
     """A Wire in build123d is a topological entity representing a connected sequence
@@ -3733,6 +3836,40 @@ class Wire(Mixin1D[TopoDS_Wire]):
 
         return ordered_edges
 
+    def geom_equal(
+        self,
+        other: Wire,
+        tol: float = 1e-6,
+        num_interpolation_points: int = 5,
+    ) -> bool:
+        """Compare two wires for geometric equality within tolerance.
+
+        This compares the geometric properties of two wires by comparing their
+        constituent edges pairwise. Two independently created wires with the
+        same geometry will return True.
+
+        Args:
+            other: Wire to compare with
+            tol: Tolerance for numeric comparisons. Defaults to 1e-6.
+            num_interpolation_points: Number of points to sample for unknown
+                curve types. Defaults to 5.
+
+        Returns:
+            bool: True if wires are geometrically equal within tolerance
+        """
+        if not isinstance(other, Wire):
+            return False
+
+        # Use order_edges to ensure consistent edge ordering and orientation
+        edges1 = self.order_edges()
+        edges2 = other.order_edges()
+        if len(edges1) != len(edges2):
+            return False
+        return all(
+            e1.geom_equal(e2, tol, num_interpolation_points)
+            for e1, e2 in zip(edges1, edges2)
+        )
+
     def param_at_point(self, point: VectorLike) -> float:
         """
         Return the normalized wire parameter for the point closest to this wire.

src/build123d/topology/shape_core.py

@@ -85,7 +85,6 @@ from OCP.BRepAlgoAPI import (
 from OCP.BRepBuilderAPI import (
     BRepBuilderAPI_Copy,
     BRepBuilderAPI_GTransform,
-    BRepBuilderAPI_MakeEdge,
     BRepBuilderAPI_MakeFace,
     BRepBuilderAPI_MakeVertex,
     BRepBuilderAPI_RightCorner,
@@ -103,7 +102,6 @@ from OCP.BRepMesh import BRepMesh_IncrementalMesh
 from OCP.BRepPrimAPI import BRepPrimAPI_MakeHalfSpace
 from OCP.BRepTools import BRepTools, BRepTools_WireExplorer
 from OCP.gce import gce_MakeLin
-from OCP.Geom import Geom_Line
 from OCP.GeomAPI import GeomAPI_ProjectPointOnSurf
 from OCP.GeomLib import GeomLib_IsPlanarSurface
 from OCP.gp import gp_Ax1, gp_Ax2, gp_Ax3, gp_Dir, gp_Pnt, gp_Trsf, gp_Vec, gp_XYZ
@@ -1341,67 +1339,89 @@ class Shape(NodeMixin, Generic[TOPODS]):
         )
 
     def intersect(
-        self, *to_intersect: Shape | Vector | Location | Axis | Plane
-    ) -> None | ShapeList[Self]:
-        """Intersection of the arguments and this shape
+        self,
+        *to_intersect: Shape | Vector | Location | Axis | Plane,
+        tolerance: float = 1e-6,
+        include_touched: bool = False,
+    ) -> ShapeList | None:
+        """Find where bodies/interiors meet (overlap or crossing geometry).
+
+        This is the main entry point for intersection operations. Handles
+        geometry conversion and delegates to subclass _intersect() implementations.
+
+        Semantics:
+            - Multiple arguments use AND (chaining): c.intersect(s1, s2) = c ∩ s1 ∩ s2
+            - Compound arguments use OR (distribution): c.intersect(Compound([s1, s2]))
+              = (c ∩ s1) ∪ (c ∩ s2)
 
         Args:
-            to_intersect (sequence of Union[Shape, Axis, Plane]): Shape(s) to
-                intersect with
+            to_intersect: Shape(s) or geometry objects to intersect with
+            tolerance: tolerance for intersection detection
+            include_touched: if True, include boundary contacts without interior
+                overlap (only relevant when Solids are involved)
 
         Returns:
-            None | ShapeList[Self]: Resulting ShapeList may contain different class
-                than self
+            ShapeList of intersection results, or None if no intersection
         """
 
-        def _to_vertex(vec: Vector) -> Vertex:
-            """Helper method to convert vector to shape"""
-            return self.__class__.cast(
-                downcast(
-                    BRepBuilderAPI_MakeVertex(gp_Pnt(vec.X, vec.Y, vec.Z)).Vertex()
-                )
-            )
+        if not to_intersect:
+            return None
 
-        def _to_edge(axis: Axis) -> Edge:
-            """Helper method to convert axis to shape"""
-            return self.__class__.cast(
-                BRepBuilderAPI_MakeEdge(
-                    Geom_Line(
-                        axis.position.to_pnt(),
-                        axis.direction.to_dir(),
-                    )
-                ).Edge()
-            )
-
-        def _to_face(plane: Plane) -> Face:
-            """Helper method to convert plane to shape"""
-            return self.__class__.cast(BRepBuilderAPI_MakeFace(plane.wrapped).Face())
-
-        # Convert any geometry objects into their respective topology objects
-        objs = []
+        # Validate input types
         for obj in to_intersect:
-            if isinstance(obj, Vector):
-                objs.append(_to_vertex(obj))
-            elif isinstance(obj, Axis):
-                objs.append(_to_edge(obj))
-            elif isinstance(obj, Plane):
-                objs.append(_to_face(obj))
-            elif isinstance(obj, Location):
-                if obj.wrapped is None:
-                    raise ValueError("Cannot intersect with an empty location")
-                objs.append(_to_vertex(tcast(Vector, obj.position)))
-            else:
-                objs.append(obj)
+            if not isinstance(obj, (Shape, Vector, Location, Axis, Plane)):
+                raise ValueError(f"Unsupported type for intersect: {type(obj)}")
 
-        # Find the shape intersections
-        intersect_op = BRepAlgoAPI_Common()
-        intersections = self._bool_op((self,), objs, intersect_op)
-        if isinstance(intersections, ShapeList):
-            return intersections or None
-        if isinstance(intersections, Shape) and not intersections.is_null:
-            return ShapeList([intersections])
+        # Chained iteration for AND semantics: c.intersect(s1, s2) = c ∩ s1 ∩ s2
+        # Geometry objects (Vector, Location, Axis, Plane) are converted in _intersect
+        common_set = ShapeList([self])
+        for other in to_intersect:
+            next_set: ShapeList = ShapeList()
+            for obj in common_set:
+                result = obj._intersect(other, tolerance, include_touched)
+                if result:
+                    next_set.extend(result.expand())
+            if not next_set:
+                return None  # AND semantics: if any step fails, no intersection
+            common_set = ShapeList(set(next_set))  # deduplicate
+        return common_set if common_set else None
+
+    def _intersect(
+        self,
+        other: Shape | Vector | Location | Axis | Plane,
+        tolerance: float = 1e-6,
+        include_touched: bool = False,
+    ) -> ShapeList | None:
+        """Single-object intersection implementation.
+
+        Base implementation returns None. Subclasses (Vertex, Mixin1D, Mixin2D,
+        Mixin3D, Compound) override this to provide actual intersection logic.
+
+        Args:
+            other: Shape or geometry object to intersect with
+            tolerance: tolerance for intersection detection
+            include_touched: if True, include boundary contacts
+
+        Returns:
+            ShapeList of intersection shapes, or None if no intersection
+        """
         return None
 
+    def touch(self, other: Shape, tolerance: float = 1e-6) -> ShapeList:
+        """Find boundary contacts between this shape and another.
+
+        Base implementation returns empty ShapeList. Subclasses (Mixin2D, Mixin3D,
+        Compound) override this to provide actual touch detection.
+
+        Args:
+            other: Shape to find contacts with
+            tolerance: tolerance for contact detection
+
+        Returns:
+            ShapeList of contact shapes (empty for base implementation)
+        """
+        return ShapeList()
+
     def is_equal(self, other: Shape) -> bool:
         """Returns True if two shapes are equal, i.e. if they share the same
         TShape with the same Locations and Orientations. Also see
@@ -2269,7 +2289,7 @@ class Shape(NodeMixin, Generic[TOPODS]):
         args: Iterable[Shape],
         tools: Iterable[Shape],
         operation: BRepAlgoAPI_BooleanOperation | BRepAlgoAPI_Splitter,
-    ) -> Self | ShapeList[Self]:
+    ) -> Self | ShapeList:
         """Generic boolean operation
 
         Args:
@@ -2279,6 +2299,7 @@ class Shape(NodeMixin, Generic[TOPODS]):
           BRepAlgoAPI_Splitter]:
 
         Returns:
+            Shape or ShapeList depending on result
 
         """
         args = list(args)
@@ -2340,6 +2361,33 @@ class Shape(NodeMixin, Generic[TOPODS]):
 
         return result
 
+    def _bool_op_list(
+        self,
+        args: Iterable[Shape],
+        tools: Iterable[Shape],
+        operation: BRepAlgoAPI_BooleanOperation | BRepAlgoAPI_Splitter,
+    ) -> ShapeList:
+        """Generic boolean operation that always returns ShapeList.
+
+        Wrapper around _bool_op that guarantees ShapeList return type,
+        wrapping single results and returning empty ShapeList for null results.
+
+        Args:
+          args: Iterable[Shape]:
+          tools: Iterable[Shape]:
+          operation: Union[BRepAlgoAPI_BooleanOperation, BRepAlgoAPI_Splitter]:
+
+        Returns:
+            ShapeList (possibly empty)
+
+        """
+        result = self._bool_op(args, tools, operation)
+        if isinstance(result, ShapeList):
+            return result
+        if result.is_null:
+            return ShapeList()
+        return ShapeList([result])
+
     def _ocp_section(
         self: Shape, other: Vertex | Edge | Wire | Face
     ) -> tuple[ShapeList[Vertex], ShapeList[Edge]]:
@@ -2607,6 +2655,29 @@ class ShapeList(list[T]):
         """Differences between two ShapeLists operator -"""
         return ShapeList(set(self) - set(other))
 
+    def expand(self) -> ShapeList:
+        """Expand by dissolving compounds, wires, and shells, filtering nulls.
+
+        Returns:
+            ShapeList with compounds dissolved to children, wires to edges,
+            shells to faces, and nulls filtered out
+        """
+        expanded: ShapeList = ShapeList()
+        for shape in self:
+            if isinstance(shape, Vector):
+                expanded.append(shape)
+            elif hasattr(shape, "wrapped"):
+                if isinstance(shape.wrapped, TopoDS_Compound):
+                    # Recursively expand nested compounds
+                    expanded.extend(ShapeList(list(shape)).expand())
+                elif isinstance(shape.wrapped, TopoDS_Shell):
+                    expanded.extend(shape.faces())
+                elif isinstance(shape.wrapped, TopoDS_Wire):
+                    expanded.extend(shape.edges())
+                elif not shape.is_null:
+                    expanded.append(shape)
+        return expanded
+
     def center(self) -> Vector:
         """The average of the center of objects within the ShapeList"""
         if not self:

src/build123d/topology/three_d.py

@@ -62,6 +62,7 @@ from typing_extensions import Self
 import OCP.TopAbs as ta
 from OCP.BRepAlgoAPI import BRepAlgoAPI_Common, BRepAlgoAPI_Cut, BRepAlgoAPI_Section
 from OCP.BRepBuilderAPI import BRepBuilderAPI_MakeSolid
+from OCP.BRepExtrema import BRepExtrema_DistShapeShape
 from OCP.BRepClass3d import BRepClass3d_SolidClassifier
 from OCP.BRepFeat import BRepFeat_MakeDPrism
 from OCP.BRepFilletAPI import BRepFilletAPI_MakeChamfer, BRepFilletAPI_MakeFillet
@@ -96,7 +97,8 @@ from OCP.TopoDS import (
     TopoDS_Solid,
     TopoDS_Wire,
 )
-from OCP.gp import gp_Ax2, gp_Pnt
+from OCP.gp import gp_Ax2, gp_Pnt, gp_Vec
+from OCP.BRepGProp import BRepGProp_Face
 from build123d.build_enums import CenterOf, GeomType, Keep, Kind, Transition, Until
 from build123d.geometry import (
     DEG2RAD,
@@ -425,131 +427,92 @@ class Mixin3D(Shape[TOPODS]):
 
         return return_value
 
-    def intersect(
-        self, *to_intersect: Shape | Vector | Location | Axis | Plane
-    ) -> None | ShapeList[Vertex | Edge | Face | Solid]:
-        """Intersect Solid with Shape or geometry object
+    def _intersect(
+        self,
+        other: Shape | Vector | Location | Axis | Plane,
+        tolerance: float = 1e-6,
+        include_touched: bool = False,
+    ) -> ShapeList | None:
+        """Single-object intersection for Solid.
+
+        Returns same-dimension overlap or crossing geometry:
+        - Solid + Solid → Solid (volume overlap)
+        - Solid + Face → Face (portion in/on solid)
+        - Solid + Edge → Edge (portion through solid)
 
         Args:
-            to_intersect (Shape | Vector | Location | Axis | Plane): objects to intersect
-
-        Returns:
-            ShapeList[Vertex | Edge | Face | Solid] | None: ShapeList of vertices, edges,
-                faces, and/or solids.
+            other: Shape or geometry object to intersect with
+            tolerance: tolerance for intersection detection
+            include_touched: if True, include boundary contacts
+                (shapes touching the solid's surface without penetrating)
         """
+        # Convert geometry objects to shapes
+        if isinstance(other, Vector):
+            other = Vertex(other)
+        elif isinstance(other, Location):
+            other = Vertex(other.position)
+        elif isinstance(other, Axis):
+            other = Edge(other)
+        elif isinstance(other, Plane):
+            other = Face(other)
 
-        def to_vector(objs: Iterable) -> ShapeList:
-            return ShapeList([Vector(v) if isinstance(v, Vertex) else v for v in objs])
-
-        def to_vertex(objs: Iterable) -> ShapeList:
-            return ShapeList([Vertex(v) if isinstance(v, Vector) else v for v in objs])
-
-        def bool_op(
-            args: Sequence,
-            tools: Sequence,
-            operation: BRepAlgoAPI_Common | BRepAlgoAPI_Section,
-        ) -> ShapeList:
-            # Wrap Shape._bool_op for corrected output
-            intersections: Shape | ShapeList = Shape()._bool_op(args, tools, operation)
-            if isinstance(intersections, ShapeList):
-                return intersections or ShapeList()
-            if isinstance(intersections, Shape) and not intersections.is_null:
-                return ShapeList([intersections])
-            return ShapeList()
-
-        def filter_shapes_by_order(shapes: ShapeList, orders: list) -> ShapeList:
-            # Remove lower order shapes from list which *appear* to be part of
-            # a higher order shape using a lazy distance check
-            # (sufficient for vertices, may be an issue for higher orders)
-            order_groups = []
-            for order in orders:
-                order_groups.append(
-                    ShapeList([s for s in shapes if isinstance(s, order)])
-                )
-
-            filtered_shapes = order_groups[-1]
-            for i in range(len(order_groups) - 1):
-                los = order_groups[i]
-                his: list = sum(order_groups[i + 1 :], [])
-                filtered_shapes.extend(
-                    ShapeList(
-                        lo
-                        for lo in los
-                        if all(lo.distance_to(hi) > TOLERANCE for hi in his)
-                    )
-                )
-
-            return filtered_shapes
-
-        common_set: ShapeList[Vertex | Edge | Face | Solid] = ShapeList([self])
-        target: Shape
-        for other in to_intersect:
-            # Conform target type
-            match other:
-                case Axis():
-                    # BRepAlgoAPI_Section seems happier if Edge isnt infinite
-                    bbox = self.bounding_box()
-                    dist = self.distance_to(other.position)
-                    dist = dist if dist >= 1 else 1
-                    target = Edge.make_line(
-                        other.position - other.direction * bbox.diagonal * dist,
-                        other.position + other.direction * bbox.diagonal * dist,
-                    )
-                case Plane():
-                    target = Face(other)
-                case Vector():
-                    target = Vertex(other)
-                case Location():
-                    target = Vertex(other.position)
-                case _ if issubclass(type(other), Shape):
-                    target = other
-                case _:
-                    raise ValueError(f"Unsupported type to_intersect: {type(other)}")
-
-            # Find common matches
-            common: list[Vertex | Edge | Wire | Face | Shell | Solid] = []
-            result: ShapeList | None
-            for obj in common_set:
-                match (obj, target):
-                    case (_, Vertex() | Edge() | Wire() | Face() | Shell() | Solid()):
-                        operation: BRepAlgoAPI_Section | BRepAlgoAPI_Common = (
-                            BRepAlgoAPI_Section()
+        def filter_redundant_touches(items: ShapeList) -> ShapeList:
+            """Remove vertices/edges that lie on higher-dimensional results."""
+            edges = [r for r in items if isinstance(r, Edge)]
+            faces = [r for r in items if isinstance(r, Face)]
+            solids = [r for r in items if isinstance(r, Solid)]
+            return ShapeList(
+                r
+                for r in items
+                if not (
+                    isinstance(r, Vertex)
+                    and (
+                        any(e.distance_to(r) <= tolerance for e in edges)
+                        or any(f.distance_to(r) <= tolerance for f in faces)
+                        or any(
+                            sf.distance_to(r) <= tolerance
+                            for s in solids
+                            for sf in s.faces()
                         )
-                        result = bool_op((obj,), (target,), operation)
-                        if (
-                            not isinstance(obj, Edge | Wire)
-                            and not isinstance(target, (Edge | Wire))
-                        ) or (isinstance(obj, Solid) or isinstance(target, Solid)):
-                            # Face + Edge combinations may produce an intersection
-                            # with Common but always with Section.
-                            # No easy way to deduplicate
-                            # Many Solid + Edge combinations need Common
-                            operation = BRepAlgoAPI_Common()
-                            result.extend(bool_op((obj,), (target,), operation))
-
-                    case _ if issubclass(type(target), Shape):
-                        result = target.intersect(obj)
-
-                if result:
-                    common.extend(result)
-
-            if common:
-                common_set = ShapeList()
-                for shape in common:
-                    if isinstance(shape, Wire):
-                        common_set.extend(shape.edges())
-                    elif isinstance(shape, Shell):
-                        common_set.extend(shape.faces())
-                    else:
-                        common_set.append(shape)
-                common_set = to_vertex(set(to_vector(common_set)))
-                common_set = filter_shapes_by_order(
-                    common_set, [Vertex, Edge, Face, Solid]
+                    )
                 )
-            else:
-                return None
+                and not (
+                    isinstance(r, Edge)
+                    and any(f.distance_to(r.center()) <= tolerance for f in faces)
+                )
+            )
 
-        return ShapeList(common_set)
+        results: ShapeList = ShapeList()
+
+        # Trim infinite edges before OCCT operations
+        if isinstance(other, Edge) and other.is_infinite:
+            bbox = self.bounding_box(optimal=False)
+            other = other.trim_infinite(
+                bbox.diagonal + (other.center() - bbox.center()).length
+            )
+
+        # Solid + Solid/Face/Shell/Edge/Wire: use Common
+        if isinstance(other, (Solid, Face, Shell, Edge, Wire)):
+            intersection = self._bool_op_list((self,), (other,), BRepAlgoAPI_Common())
+            results.extend(intersection.expand())
+        # Solid + Vertex: point containment check
+        elif isinstance(other, Vertex):
+            if self.is_inside(Vector(other), tolerance):
+                results.append(other)
+
+        # Delegate to higher-order shapes (Compound)
+        # Don't pass include_touched - outer caller handles touches
+        else:
+            result = other._intersect(self, tolerance, include_touched=False)
+            if result:
+                results.extend(result)
+
+        # Add boundary contacts if requested (only Solid has touch method)
+        if include_touched and isinstance(self, Solid):
+            results.extend(self.touch(other, tolerance))
+            results = filter_redundant_touches(ShapeList(set(results)))
+
+        return results if results else None
 
     def is_inside(self, point: VectorLike, tolerance: float = 1.0e-6) -> bool:
         """Returns whether or not the point is inside a solid or compound
@@ -794,6 +757,210 @@ class Solid(Mixin3D[TopoDS_Solid]):
         # when density == 1, mass == volume
         return Shape.compute_mass(self)
 
+    # ---- Instance Methods ----
+
+    def touch(
+        self,
+        other: Shape,
+        tolerance: float = 1e-6,
+        found_solids: ShapeList | None = None,
+    ) -> ShapeList[Vertex | Edge | Face]:
+        """Find where this Solid's boundary contacts another shape.
+
+        Returns geometry where boundaries contact without interior overlap:
+        - Solid + Solid → Face + Edge + Vertex (all boundary contacts)
+        - Solid + Face/Shell → Face + Edge + Vertex (boundary contacts)
+        - Solid + Edge/Wire → Vertex (edge endpoints on solid boundary)
+        - Solid + Vertex → Vertex if on boundary
+        - Solid + Compound → distributes over compound elements
+
+        Args:
+            other: Shape to check boundary contacts with
+            tolerance: tolerance for contact detection
+            found_solids: pre-found intersection solids to filter against
+
+        Returns:
+            ShapeList of boundary contact geometry (empty if no contact)
+        """
+
+        # Helper functions for common geometric checks (for readability)
+        # Single shape versions for checking against one shapes
+        def vertex_on_edge(v: Vertex, e: Edge) -> bool:
+            return v.distance_to(e) <= tolerance
+
+        def vertex_on_face(v: Vertex, f: Face) -> bool:
+            return v.distance_to(f) <= tolerance
+
+        def edge_on_face(e: Edge, f: Face) -> bool:
+            # Can't use distance_to (e.g. normal vector would match), need Common
+            return bool(self._bool_op_list((e,), (f,), BRepAlgoAPI_Common()))
+
+        # Multi shape versions for checking against multiple shapes
+        def vertex_on_edges(v: Vertex, edges: Iterable[Edge]) -> bool:
+            return any(vertex_on_edge(v, e) for e in edges)
+
+        def vertex_on_faces(v: Vertex, faces: Iterable[Face]) -> bool:
+            return any(vertex_on_face(v, f) for f in faces)
+
+        def edge_on_faces(e: Edge, faces: Iterable[Face]) -> bool:
+            return any(edge_on_face(e, f) for f in faces)
+
+        def face_point_normal(face: Face, u: float, v: float) -> tuple[Vector, Vector]:
+            """Get both position and normal at UV coordinates.
+            Args
+                u (float): the horizontal coordinate in the parameter space of the Face,
+                    between 0.0 and 1.0
+                v (float): the vertical coordinate in the parameter space of the Face,
+                    between 0.0 and 1.0
+            Returns:
+                tuple[Vector, Vector]: [point on Face, normal at point]
+            """
+            u0, u1, v0, v1 = face._uv_bounds()
+            u_val = u0 + u * (u1 - u0)
+            v_val = v0 + v * (v1 - v0)
+            gp_pnt = gp_Pnt()
+            gp_norm = gp_Vec()
+            BRepGProp_Face(face.wrapped).Normal(u_val, v_val, gp_pnt, gp_norm)
+            return Vector(gp_pnt), Vector(gp_norm)
+
+        def faces_equal(f1: Face, f2: Face, grid_size: int = 4) -> bool:
+            """Check if two faces are geometrically equal.
+
+            Face == uses topological equality (same OCC object), but we need
+            geometric equality. For performance reasons apply a heuristic
+            approach: Compare a grid of UV sample points, checking both position and
+            normal direction match within tolerance.
+            """
+            # Early reject: bounding box check
+            bb1 = f1.bounding_box(optimal=False)
+            bb2 = f2.bounding_box(optimal=False)
+            if not bb1.overlaps(bb2, tolerance):
+                return False
+
+            # Compare grid_size x grid_size grid of points in UV space
+            for i in range(grid_size):
+                u = i / (grid_size - 1)
+                for j in range(grid_size):
+                    v = j / (grid_size - 1)
+                    pos1, norm1 = face_point_normal(f1, u, v)
+                    pos2, norm2 = face_point_normal(f2, u, v)
+                    if (pos1 - pos2).length > tolerance or abs(norm1.dot(norm2)) < 0.99:
+                        return False
+            return True
+
+        def is_duplicate(shape: Shape, existing: Iterable[Shape]) -> bool:
+            if isinstance(shape, Vertex):
+                return any(
+                    isinstance(v, Vertex) and Vector(shape) == Vector(v)
+                    for v in existing
+                )
+            if isinstance(shape, Edge):
+                return any(
+                    isinstance(e, Edge) and shape.geom_equal(e, tolerance)
+                    for e in existing
+                )
+            if isinstance(shape, Face):
+                # Heuristic approach
+                return any(
+                    isinstance(f, Face) and faces_equal(shape, f) for f in existing
+                )
+            return False
+
+        results: ShapeList = ShapeList()
+
+        if isinstance(other, (Solid, Face, Shell)):
+            # Unified handling: iterate over face pairs
+            # For Solid+Solid: get intersection solids to filter results that bound them
+            intersect_faces = []
+            if isinstance(other, Solid):
+                if found_solids is None:
+                    found_solids = ShapeList(
+                        self._intersect(other, tolerance, include_touched=False) or []
+                    )
+                intersect_faces = [f for s in found_solids for f in s.faces()]
+
+            # Pre-calculate bounding boxes for early rejection
+            self_faces = [(f, f.bounding_box(optimal=False)) for f in self.faces()]
+            other_faces = [(f, f.bounding_box(optimal=False)) for f in other.faces()]
+
+            # First pass: collect touch/intersect results from face pairs,
+            # filtering against intersection solid faces
+            raw_results: ShapeList = ShapeList()
+            for sf, sf_bb in self_faces:
+                for of, of_bb in other_faces:
+                    if not sf_bb.overlaps(of_bb, tolerance):
+                        continue
+
+                    # Process touch first (cheap), then intersect (expensive)
+                    # Face touch gives tangent vertices
+                    for r in sf.touch(of, tolerance=tolerance):
+                        if not is_duplicate(r, raw_results) and not vertex_on_faces(
+                            r, intersect_faces
+                        ):
+                            raw_results.append(r)
+
+                    # Face intersect gives shared faces/edges (touch handled above)
+                    for r in sf.intersect(of, tolerance=tolerance) or []:
+                        if not is_duplicate(r, raw_results) and not edge_on_faces(
+                            r, intersect_faces
+                        ):
+                            raw_results.append(r)
+
+            # Second pass: filter lower-dimensional results against higher-dimensional
+            all_faces = [f for f in raw_results if isinstance(f, Face)]
+            all_edges = [e for e in raw_results if isinstance(e, Edge)]
+            for r in raw_results:
+                if (
+                    isinstance(r, Face)
+                    or (isinstance(r, Edge) and not edge_on_faces(r, all_faces))
+                    or (
+                        isinstance(r, Vertex)
+                        and not vertex_on_faces(r, all_faces)
+                        and not vertex_on_edges(r, all_edges)
+                    )
+                ):
+                    results.append(r)
+
+        elif isinstance(other, (Edge, Wire)):
+            # Solid + Edge: find where edge endpoints touch solid boundary
+            # Pre-calculate bounding boxes (optimal=False for speed, used for filtering)
+            self_faces = [(f, f.bounding_box(optimal=False)) for f in self.faces()]
+            other_bb = other.bounding_box(optimal=False)
+
+            for ov in other.vertices():
+                for sf, _ in self_faces:
+                    if vertex_on_face(ov, sf):
+                        results.append(ov)
+                        break
+
+            # Use BRepExtrema to find all tangent contacts (edge tangent to surface)
+            for sf, sf_bb in self_faces:
+                if not sf_bb.overlaps(other_bb, tolerance):
+                    continue
+                extrema = BRepExtrema_DistShapeShape(sf.wrapped, other.wrapped)
+                if extrema.IsDone() and extrema.Value() <= tolerance:
+                    for i in range(1, extrema.NbSolution() + 1):
+                        pnt1 = extrema.PointOnShape1(i)
+                        pnt2 = extrema.PointOnShape2(i)
+                        if pnt1.Distance(pnt2) <= tolerance:
+                            new_vertex = Vertex(pnt1.X(), pnt1.Y(), pnt1.Z())
+                            if not is_duplicate(new_vertex, results):
+                                results.append(new_vertex)
+
+        elif isinstance(other, Vertex):
+            # Solid + Vertex: check if vertex is on boundary
+            for sf in self.faces():
+                if vertex_on_face(other, sf):
+                    results.append(other)
+                    break
+
+        # Delegate to other shapes (Compound iterates, others return empty)
+        else:
+            results.extend(other.touch(self, tolerance))
+
+        # Remove duplicates using Shape's __hash__ and __eq__
+        return ShapeList(set(results))
+
     # ---- Class Methods ----
 
     @classmethod

src/build123d/topology/two_d.py

@@ -69,6 +69,7 @@ from OCP.BRep import BRep_Builder, BRep_Tool
 from OCP.BRepAdaptor import BRepAdaptor_Curve, BRepAdaptor_Surface
 from OCP.BRepAlgo import BRepAlgo
 from OCP.BRepAlgoAPI import BRepAlgoAPI_Common, BRepAlgoAPI_Section
+from OCP.BRepExtrema import BRepExtrema_DistShapeShape
 from OCP.BRepBuilderAPI import (
     BRepBuilderAPI_MakeEdge,
     BRepBuilderAPI_MakeFace,
@@ -118,7 +119,13 @@ from OCP.TColStd import (
     TColStd_HArray2OfReal,
 )
 from OCP.TopExp import TopExp
-from OCP.TopoDS import TopoDS, TopoDS_Face, TopoDS_Shape, TopoDS_Shell, TopoDS_Solid
+from OCP.TopoDS import (
+    TopoDS,
+    TopoDS_Face,
+    TopoDS_Shape,
+    TopoDS_Shell,
+    TopoDS_Solid,
+)
 from OCP.TopTools import TopTools_IndexedDataMapOfShapeListOfShape, TopTools_ListOfShape
 from ocp_gordon import interpolate_curve_network
 from typing_extensions import Self
@@ -277,127 +284,217 @@ class Mixin2D(ABC, Shape[TOPODS]):
 
         return result
 
-    def intersect(
-        self, *to_intersect: Shape | Vector | Location | Axis | Plane
-    ) -> None | ShapeList[Vertex | Edge | Face]:
-        """Intersect Face with Shape or geometry object
+    def _intersect(
+        self,
+        other: Shape | Vector | Location | Axis | Plane,
+        tolerance: float = 1e-6,
+        include_touched: bool = False,
+    ) -> ShapeList | None:
+        """Single-object intersection for Face/Shell.
+
+        Returns same-dimension overlap or crossing geometry:
+        - 2D + 2D → Face (coplanar overlap) + Edge (crossing curves)
+        - 2D + Edge → Edge (on surface) + Vertex (piercing)
+        - 2D + Solid/Compound → delegates to other._intersect(self)
 
         Args:
-            to_intersect (Shape | Vector | Location | Axis | Plane): objects to intersect
+            other: Shape or geometry object to intersect with
+            tolerance: tolerance for intersection detection
+            include_touched: if True, include boundary contacts
+                (only relevant when Solids are involved)
+        """
+        # Convert geometry objects to shapes
+        if isinstance(other, Vector):
+            other = Vertex(other)
+        elif isinstance(other, Location):
+            other = Vertex(other.position)
+        elif isinstance(other, Axis):
+            other = Edge(other)
+        elif isinstance(other, Plane):
+            other = Face(other)
+
+        def filter_edges(
+            section_edges: ShapeList[Edge], common_edges: ShapeList[Edge]
+        ) -> ShapeList[Edge]:
+            """Filter section edges, keeping only edges not on common face boundaries."""
+            # Pre-compute bounding boxes for both sets (optimal=False for speed, filtering only)
+            section_bboxes = [(e, e.bounding_box(optimal=False)) for e in section_edges]
+            common_bboxes = [
+                (ce, ce.bounding_box(optimal=False)) for ce in common_edges
+            ]
+
+            # Filter: remove section edges that coincide with common face boundaries
+            filtered: ShapeList = ShapeList()
+            for edge, edge_bbox in section_bboxes:
+                is_common = any(
+                    edge_bbox.overlaps(ce_bbox, tolerance)
+                    and edge.distance_to(ce) <= tolerance
+                    for ce, ce_bbox in common_bboxes
+                )
+                if not is_common:
+                    filtered.append(edge)
+            return filtered
+
+        results: ShapeList = ShapeList()
+
+        # Trim infinite edges before OCCT operations
+        if isinstance(other, Edge) and other.is_infinite:
+            bbox = self.bounding_box(optimal=False)
+            other = other.trim_infinite(
+                bbox.diagonal + (other.center() - bbox.center()).length
+            )
+
+        # 2D + 2D: Common (coplanar overlap) AND Section (crossing curves)
+        if isinstance(other, (Face, Shell)):
+            # Common for coplanar overlap
+            common = self._bool_op_list((self,), (other,), BRepAlgoAPI_Common())
+            common_faces = common.expand()
+            results.extend(common_faces)
+
+            # Section for crossing curves (only edges, not vertices)
+            # Vertices from Section are boundary contacts (touch), not intersections
+            section = self._bool_op_list((self,), (other,), BRepAlgoAPI_Section())
+            section_edges = ShapeList(
+                [s for s in section if isinstance(s, Edge)]
+            ).expand()
+
+            if not common_faces:
+                # No coplanar overlap - all section edges are valid crossings
+                results.extend(section_edges)
+            else:
+                # Filter out edges on common face boundaries
+                # (Section returns boundary of overlap region which are not crossings)
+                common_edges: ShapeList[Edge] = ShapeList()
+                for face in common_faces:
+                    common_edges.extend(face.edges())
+                results.extend(filter_edges(section_edges, common_edges))
+
+        # 2D + Edge: Section for intersection
+        elif isinstance(other, (Edge, Wire)):
+            section = self._bool_op_list((self,), (other,), BRepAlgoAPI_Section())
+            results.extend(section)
+
+        # 2D + Vertex: point containment on surface
+        elif isinstance(other, Vertex):
+            if other.distance_to(self) <= tolerance:
+                results.append(other)
+
+        # Delegate to higher-order shapes (Solid, etc.)
+        else:
+            result = other._intersect(self, tolerance, include_touched)
+            if result:
+                results.extend(result)
+
+        # Add boundary contacts if requested
+        if include_touched and isinstance(other, (Face, Shell)):
+            found_faces = ShapeList(r for r in results if isinstance(r, Face))
+            found_edges = ShapeList(r for r in results if isinstance(r, Edge))
+            results.extend(self.touch(other, tolerance, found_faces, found_edges))
+
+        return results if results else None
+
+    def touch(
+        self,
+        other: Shape,
+        tolerance: float = 1e-6,
+        found_faces: ShapeList | None = None,
+        found_edges: ShapeList | None = None,
+    ) -> ShapeList:
+        """Find boundary contacts between this 2D shape and another shape.
+
+        Returns the highest-dimensional contact at each location, filtered to
+        avoid returning lower-dimensional boundaries of higher-dimensional contacts.
+
+        For Face/Shell:
+        - Face + Face → Vertex (shared corner or crossing point without edge/face overlap)
+        - Face + Edge/Vertex → no touch (intersect already returns dim 0)
+
+        Args:
+            other: Shape to find contacts with
+            tolerance: tolerance for contact detection
+            found_faces: pre-found faces to filter against (from Mixin3D.touch)
+            found_edges: pre-found edges to filter against (from Mixin3D.touch)
 
         Returns:
-            ShapeList[Vertex | Edge | Face] | None: ShapeList of vertices, edges, and/or
-                faces.
+            ShapeList of contact shapes (Vertex only for 2D+2D)
         """
 
-        def to_vector(objs: Iterable) -> ShapeList:
-            return ShapeList([Vector(v) if isinstance(v, Vertex) else v for v in objs])
+        # Helper functions for common geometric checks
+        def vertex_on_edges(v: Vertex, edges: Iterable[Edge]) -> bool:
+            return any(v.distance_to(e) <= tolerance for e in edges)
 
-        def to_vertex(objs: Iterable) -> ShapeList:
-            return ShapeList([Vertex(v) if isinstance(v, Vector) else v for v in objs])
+        def vertex_on_faces(v: Vertex, faces: Iterable[Face]) -> bool:
+            return any(v.distance_to(f) <= tolerance for f in faces)
 
-        def bool_op(
-            args: Sequence,
-            tools: Sequence,
-            operation: BRepAlgoAPI_Section | BRepAlgoAPI_Common,
-        ) -> ShapeList:
-            # Wrap Shape._bool_op for corrected output
-            intersections: Shape | ShapeList = Shape()._bool_op(args, tools, operation)
-            if isinstance(intersections, ShapeList):
-                return intersections or ShapeList()
-            if isinstance(intersections, Shape) and not intersections.is_null:
-                return ShapeList([intersections])
-            return ShapeList()
+        def is_duplicate(v: Vertex, vertices: Iterable[Vertex]) -> bool:
+            vec = Vector(v)
+            return any(vec == Vector(ov) for ov in vertices)
 
-        def filter_shapes_by_order(shapes: ShapeList, orders: list) -> ShapeList:
-            # Remove lower order shapes from list which *appear* to be part of
-            # a higher order shape using a lazy distance check
-            # (sufficient for vertices, may be an issue for higher orders)
-            order_groups = []
-            for order in orders:
-                order_groups.append(
-                    ShapeList([s for s in shapes if isinstance(s, order)])
+        results: ShapeList = ShapeList()
+
+        if isinstance(other, (Face, Shell)):
+            # Get intersect results to filter against if not provided (direct call)
+            if found_faces is None:
+                found_faces = ShapeList()
+                found_edges = ShapeList()
+                intersect_results = self._intersect(
+                    other, tolerance, include_touched=False
                 )
+                if intersect_results:
+                    for r in intersect_results:
+                        if isinstance(r, Face):
+                            found_faces.append(r)
+                        elif isinstance(r, Edge):
+                            found_edges.append(r)
+            elif found_edges is None:  # for mypy
+                found_edges = ShapeList()
 
-            filtered_shapes = order_groups[-1]
-            for i in range(len(order_groups) - 1):
-                los = order_groups[i]
-                his: list = sum(order_groups[i + 1 :], [])
-                filtered_shapes.extend(
-                    ShapeList(
-                        lo
-                        for lo in los
-                        if all(lo.distance_to(hi) > TOLERANCE for hi in his)
-                    )
-                )
+            # Use BRepExtrema to find all contact points
+            # (vertex-vertex, vertex-edge, vertex-face)
+            found_vertices: ShapeList = ShapeList()
+            extrema = BRepExtrema_DistShapeShape()
+            extrema.SetDeflection(
+                tolerance * 1e-3
+            )  # Higher precision to avoid duplicate solutions
+            extrema.LoadS1(self.wrapped)
+            extrema.LoadS2(other.wrapped)
+            extrema.Perform()
+            if extrema.IsDone() and extrema.Value() <= tolerance:
+                for i in range(1, extrema.NbSolution() + 1):
+                    pnt1 = extrema.PointOnShape1(i)
+                    pnt2 = extrema.PointOnShape2(i)
+                    if pnt1.Distance(pnt2) > tolerance:
+                        continue
 
-            return filtered_shapes
+                    new_vertex = Vertex(pnt1.X(), pnt1.Y(), pnt1.Z())
 
-        common_set: ShapeList[Vertex | Edge | Face | Shell] = ShapeList([self])
-        target: Shape
-        for other in to_intersect:
-            # Conform target type
-            match other:
-                case Axis():
-                    # BRepAlgoAPI_Section seems happier if Edge isnt infinite
-                    bbox = self.bounding_box()
-                    dist = self.distance_to(other.position)
-                    dist = dist if dist >= 1 else 1
-                    target = Edge.make_line(
-                        other.position - other.direction * bbox.diagonal * dist,
-                        other.position + other.direction * bbox.diagonal * dist,
-                    )
-                case Plane():
-                    target = Face(other)
-                case Vector():
-                    target = Vertex(other)
-                case Location():
-                    target = Vertex(other.position)
-                case _ if issubclass(type(other), Shape):
-                    target = other
-                case _:
-                    raise ValueError(f"Unsupported type to_intersect: {type(other)}")
+                    # Skip duplicates early (cheap check)
+                    if is_duplicate(new_vertex, found_vertices):
+                        continue
 
-            # Find common matches
-            common: list[Vertex | Edge | Wire | Face | Shell] = []
-            result: ShapeList | None
-            for obj in common_set:
-                match (obj, target):
-                    case (_, Vertex() | Edge() | Wire() | Face() | Shell()):
-                        operation: BRepAlgoAPI_Section | BRepAlgoAPI_Common = (
-                            BRepAlgoAPI_Section()
-                        )
-                        result = bool_op((obj,), (target,), operation)
-                        if not isinstance(obj, Edge | Wire) and not isinstance(
-                            target, (Edge | Wire)
-                        ):
-                            # Face + Edge combinations may produce an intersection
-                            # with Common but always with Section.
-                            # No easy way to deduplicate
-                            operation = BRepAlgoAPI_Common()
-                            result.extend(bool_op((obj,), (target,), operation))
+                    # Skip edge-edge intersections, but allow corner touches
+                    if (
+                        vertex_on_edges(new_vertex, self.edges())
+                        and vertex_on_edges(new_vertex, other.edges())
+                        and not is_duplicate(new_vertex, self.vertices())
+                        and not is_duplicate(new_vertex, other.vertices())
+                    ):
+                        continue
 
-                    case _ if issubclass(type(target), Shape):
-                        result = target.intersect(obj)
+                    # Filter: only keep vertices that are not boundaries of
+                    # higher-dimensional contacts (faces or edges)
+                    if not vertex_on_faces(
+                        new_vertex, found_faces
+                    ) and not vertex_on_edges(new_vertex, found_edges):
+                        results.append(new_vertex)
+                        found_vertices.append(new_vertex)
 
-                if result:
-                    common.extend(result)
+        # Face + Edge/Vertex: no touch (intersect already covers dim 0)
+        # Delegate to other shapes (Compound iterates, others return empty)
+        else:
+            results.extend(other.touch(self, tolerance))
 
-            if common:
-                common_set = ShapeList()
-                for shape in common:
-                    if isinstance(shape, Wire):
-                        common_set.extend(shape.edges())
-                    elif isinstance(shape, Shell):
-                        common_set.extend(shape.faces())
-                    else:
-                        common_set.append(shape)
-                common_set = to_vertex(set(to_vector(common_set)))
-                common_set = filter_shapes_by_order(common_set, [Vertex, Edge, Face])
-            else:
-                return None
-
-        return ShapeList(common_set)
+        return results
 
     @abstractmethod
     def location_at(self, *args: Any, **kwargs: Any) -> Location:

src/build123d/topology/zero_d.py

@@ -168,44 +168,46 @@ class Vertex(Shape[TopoDS_Vertex]):
         """extrude - invalid operation for Vertex"""
         raise NotImplementedError("Vertices can't be created by extrusion")
 
-    def intersect(
-        self, *to_intersect: Shape | Vector | Location | Axis | Plane
-    ) -> ShapeList[Vertex] | None:
-        """Intersection of vertex and geometric objects or shapes.
+    def _intersect(
+        self,
+        other: Shape | Vector | Location | Axis | Plane,
+        tolerance: float = 1e-6,
+        include_touched: bool = False,
+    ) -> ShapeList | None:
+        """Single-object intersection for Vertex.
+
+        For a vertex (0D), intersection means the vertex lies on/in the other shape.
 
         Args:
-            to_intersect (sequence of [Shape | Vector | Location | Axis | Plane]):
-                Objects(s) to intersect with
-
-        Returns:
-            ShapeList[Vertex] | None: Vertex intersection in a ShapeList or None
+            other: Shape or geometry object to intersect with
+            tolerance: tolerance for intersection detection
+            include_touched: if True, include boundary contacts
+                (only relevant when Solids are involved)
         """
-        common = Vector(self)
-        result: Shape | ShapeList[Shape] | Vector | None
-        for obj in to_intersect:
-            # Treat as Vector, otherwise call intersection from Shape
-            match obj:
-                case Vertex():
-                    result = common.intersect(Vector(obj))
-                case Vector() | Location() | Axis() | Plane():
-                    result = obj.intersect(common)
-                case _ if issubclass(type(obj), Shape):
-                    result = obj.intersect(self)
-                case _:
-                    raise ValueError(f"Unsupported type to_intersect:: {type(obj)}")
+        # Convert geometry objects to Vertex
+        if isinstance(other, Vector):
+            other = Vertex(other)
+        elif isinstance(other, Location):
+            other = Vertex(other.position)
+        elif isinstance(other, Axis):
+            # Check if vertex lies on the axis
+            if other.intersect(self.center()):
+                return ShapeList([self])
+            return None
+        elif isinstance(other, Plane):
+            # Check if vertex lies on the plane
+            if other.contains(self.center(), tolerance):
+                return ShapeList([self])
+            return None
 
-            if isinstance(result, Vector) and result == common:
-                pass
-            elif (
-                isinstance(result, list)
-                and len(result) == 1
-                and Vector(result[0]) == common
-            ):
-                pass
-            else:
-                return None
+        if isinstance(other, Vertex):
+            # Vertex + Vertex: check distance
+            if self.distance_to(other) <= tolerance:
+                return ShapeList([self])
+            return None
 
-        return ShapeList([self])
+        # Delegate to higher-dimensional shape (including Compound)
+        return other._intersect(self, tolerance, include_touched)
 
     # ---- Instance Methods ----
 

tests/test_direct_api/test_geom_equal.py

@@ -0,0 +1,968 @@
+"""Tests for Edge.geom_equal and Wire.geom_equal methods."""
+
+import pytest
+from build123d import (
+    Vertex,
+    Edge,
+    Wire,
+    Spline,
+    Rectangle,
+    Circle,
+    Ellipse,
+    Bezier,
+    GeomType,
+    Location,
+    Plane,
+)
+
+
+class TestEdgeGeomEqualLine:
+    """Tests for Edge.geom_equal with LINE type."""
+
+    def test_same_line(self):
+        e1 = Edge.make_line((0, 0, 0), (1, 1, 1))
+        e2 = Edge.make_line((0, 0, 0), (1, 1, 1))
+        assert e1.geom_type == GeomType.LINE
+        assert e1.geom_equal(e2)
+
+    def test_different_line(self):
+        e1 = Edge.make_line((0, 0, 0), (1, 1, 1))
+        e2 = Edge.make_line((0, 0, 0), (1, 1, 2))
+        assert not e1.geom_equal(e2)
+
+
+class TestEdgeGeomEqualCircle:
+    """Tests for Edge.geom_equal with CIRCLE type."""
+
+    def test_same_circle(self):
+        c1 = Circle(10)
+        c2 = Circle(10)
+        e1 = c1.edge()
+        e2 = c2.edge()
+        assert e1.geom_type == GeomType.CIRCLE
+        assert e1.geom_equal(e2)
+
+    def test_different_radius(self):
+        c1 = Circle(10)
+        c2 = Circle(11)
+        e1 = c1.edge()
+        e2 = c2.edge()
+        assert not e1.geom_equal(e2)
+
+    def test_same_arc(self):
+        e1 = Edge.make_circle(10, start_angle=0, end_angle=90)
+        e2 = Edge.make_circle(10, start_angle=0, end_angle=90)
+        assert e1.geom_equal(e2)
+
+    def test_different_arc_angle(self):
+        e1 = Edge.make_circle(10, start_angle=0, end_angle=90)
+        e2 = Edge.make_circle(10, start_angle=0, end_angle=180)
+        assert not e1.geom_equal(e2)
+
+    def test_different_circle_from_revolve(self):
+        """Two circles with same radius/endpoints but different center/axis."""
+        from build123d import Axis, Line, RadiusArc, make_face, revolve
+
+        f1 = make_face(RadiusArc((5, 0), (-5, 0), 15) + Line((5, 0), (-5, 0)))
+        p1 = revolve(f1, Axis.X, 90)
+        value1, value2 = p1.edges().filter_by(GeomType.CIRCLE)
+        value2 = value2.reversed()
+        # These circles have same endpoints after reversal but different center/axis
+        assert not value1.geom_equal(value2)
+
+    def test_different_location(self):
+        """Circles with same radius but different center location."""
+        e1 = Edge.make_circle(10)
+        e2 = Edge.make_circle(10).locate(Location((5, 0, 0)))
+        assert not e1.geom_equal(e2)
+
+    def test_same_location(self):
+        """Circles with same radius and same non-origin location."""
+        e1 = Edge.make_circle(10).locate(Location((5, 5, 0)))
+        e2 = Edge.make_circle(10).locate(Location((5, 5, 0)))
+        assert e1.geom_equal(e2)
+
+    def test_different_axis(self):
+        """Circles with same radius but different axis direction."""
+        e1 = Edge.make_circle(10, plane=Plane.XY)
+        e2 = Edge.make_circle(10, plane=Plane.XZ)
+        assert not e1.geom_equal(e2)
+
+    def test_same_axis(self):
+        """Circles with same radius and same non-default axis."""
+        e1 = Edge.make_circle(10, plane=Plane.YZ)
+        e2 = Edge.make_circle(10, plane=Plane.YZ)
+        assert e1.geom_equal(e2)
+
+
+class TestEdgeGeomEqualEllipse:
+    """Tests for Edge.geom_equal with ELLIPSE type."""
+
+    def test_same_ellipse(self):
+        el1 = Ellipse(10, 5)
+        el2 = Ellipse(10, 5)
+        e1 = el1.edge()
+        e2 = el2.edge()
+        assert e1.geom_type == GeomType.ELLIPSE
+        assert e1.geom_equal(e2)
+
+    def test_different_major_radius(self):
+        el1 = Ellipse(10, 5)
+        el2 = Ellipse(11, 5)
+        e1 = el1.edge()
+        e2 = el2.edge()
+        assert not e1.geom_equal(e2)
+
+    def test_different_minor_radius(self):
+        el1 = Ellipse(10, 5)
+        el2 = Ellipse(10, 6)
+        e1 = el1.edge()
+        e2 = el2.edge()
+        assert not e1.geom_equal(e2)
+
+    def test_different_location(self):
+        """Ellipses with same radii but different center location."""
+        e1 = Edge.make_ellipse(10, 5)
+        e2 = Edge.make_ellipse(10, 5).locate(Location((5, 0, 0)))
+        assert not e1.geom_equal(e2)
+
+    def test_same_location(self):
+        """Ellipses with same radii and same non-origin location."""
+        e1 = Edge.make_ellipse(10, 5).locate(Location((5, 5, 0)))
+        e2 = Edge.make_ellipse(10, 5).locate(Location((5, 5, 0)))
+        assert e1.geom_equal(e2)
+
+    def test_different_axis(self):
+        """Ellipses with same radii but different axis direction."""
+        e1 = Edge.make_ellipse(10, 5, plane=Plane.XY)
+        e2 = Edge.make_ellipse(10, 5, plane=Plane.XZ)
+        assert not e1.geom_equal(e2)
+
+    def test_same_axis(self):
+        """Ellipses with same radii and same non-default axis."""
+        e1 = Edge.make_ellipse(10, 5, plane=Plane.YZ)
+        e2 = Edge.make_ellipse(10, 5, plane=Plane.YZ)
+        assert e1.geom_equal(e2)
+
+
+class TestEdgeGeomEqualHyperbola:
+    """Tests for Edge.geom_equal with HYPERBOLA type."""
+
+    def test_same_hyperbola(self):
+        e1 = Edge.make_hyperbola(10, 5, start_angle=-45, end_angle=45)
+        e2 = Edge.make_hyperbola(10, 5, start_angle=-45, end_angle=45)
+        assert e1.geom_type == GeomType.HYPERBOLA
+        assert e1.geom_equal(e2)
+
+    def test_different_x_radius(self):
+        e1 = Edge.make_hyperbola(10, 5, start_angle=-45, end_angle=45)
+        e2 = Edge.make_hyperbola(11, 5, start_angle=-45, end_angle=45)
+        assert not e1.geom_equal(e2)
+
+    def test_different_y_radius(self):
+        e1 = Edge.make_hyperbola(10, 5, start_angle=-45, end_angle=45)
+        e2 = Edge.make_hyperbola(10, 6, start_angle=-45, end_angle=45)
+        assert not e1.geom_equal(e2)
+
+    def test_different_location(self):
+        """Hyperbolas with same radii but different center location."""
+        e1 = Edge.make_hyperbola(10, 5, start_angle=-45, end_angle=45)
+        e2 = Edge.make_hyperbola(10, 5, start_angle=-45, end_angle=45).locate(
+            Location((5, 0, 0))
+        )
+        assert not e1.geom_equal(e2)
+
+    def test_same_location(self):
+        """Hyperbolas with same radii and same non-origin location."""
+        e1 = Edge.make_hyperbola(10, 5, start_angle=-45, end_angle=45).locate(
+            Location((5, 5, 0))
+        )
+        e2 = Edge.make_hyperbola(10, 5, start_angle=-45, end_angle=45).locate(
+            Location((5, 5, 0))
+        )
+        assert e1.geom_equal(e2)
+
+    def test_different_axis(self):
+        """Hyperbolas with same radii but different axis direction."""
+        e1 = Edge.make_hyperbola(10, 5, plane=Plane.XY, start_angle=-45, end_angle=45)
+        e2 = Edge.make_hyperbola(10, 5, plane=Plane.XZ, start_angle=-45, end_angle=45)
+        assert not e1.geom_equal(e2)
+
+    def test_same_axis(self):
+        """Hyperbolas with same radii and same non-default axis."""
+        e1 = Edge.make_hyperbola(10, 5, plane=Plane.YZ, start_angle=-45, end_angle=45)
+        e2 = Edge.make_hyperbola(10, 5, plane=Plane.YZ, start_angle=-45, end_angle=45)
+        assert e1.geom_equal(e2)
+
+
+class TestEdgeGeomEqualParabola:
+    """Tests for Edge.geom_equal with PARABOLA type."""
+
+    def test_same_parabola(self):
+        e1 = Edge.make_parabola(5, start_angle=0, end_angle=60)
+        e2 = Edge.make_parabola(5, start_angle=0, end_angle=60)
+        assert e1.geom_type == GeomType.PARABOLA
+        assert e1.geom_equal(e2)
+
+    def test_different_focal_length(self):
+        e1 = Edge.make_parabola(5, start_angle=0, end_angle=60)
+        e2 = Edge.make_parabola(6, start_angle=0, end_angle=60)
+        assert not e1.geom_equal(e2)
+
+    def test_different_location(self):
+        """Parabolas with same focal length but different vertex location."""
+        e1 = Edge.make_parabola(5, start_angle=0, end_angle=60)
+        e2 = Edge.make_parabola(5, start_angle=0, end_angle=60).locate(
+            Location((5, 0, 0))
+        )
+        assert not e1.geom_equal(e2)
+
+    def test_same_location(self):
+        """Parabolas with same focal length and same non-origin location."""
+        e1 = Edge.make_parabola(5, start_angle=0, end_angle=60).locate(
+            Location((5, 5, 0))
+        )
+        e2 = Edge.make_parabola(5, start_angle=0, end_angle=60).locate(
+            Location((5, 5, 0))
+        )
+        assert e1.geom_equal(e2)
+
+    def test_different_axis(self):
+        """Parabolas with same focal length but different axis direction."""
+        e1 = Edge.make_parabola(5, plane=Plane.XY, start_angle=0, end_angle=60)
+        e2 = Edge.make_parabola(5, plane=Plane.XZ, start_angle=0, end_angle=60)
+        assert not e1.geom_equal(e2)
+
+    def test_same_axis(self):
+        """Parabolas with same focal length and same non-default axis."""
+        e1 = Edge.make_parabola(5, plane=Plane.YZ, start_angle=0, end_angle=60)
+        e2 = Edge.make_parabola(5, plane=Plane.YZ, start_angle=0, end_angle=60)
+        assert e1.geom_equal(e2)
+
+
+class TestEdgeGeomEqualBezier:
+    """Tests for Edge.geom_equal with BEZIER type."""
+
+    def test_same_bezier(self):
+        pts = [(0, 0), (1, 1), (2, 0)]
+        b1 = Bezier(*pts)
+        b2 = Bezier(*pts)
+        e1 = b1.edge()
+        e2 = b2.edge()
+        assert e1.geom_type == GeomType.BEZIER
+        assert e1.geom_equal(e2)
+
+    def test_different_bezier(self):
+        b1 = Bezier((0, 0), (1, 1), (2, 0))
+        b2 = Bezier((0, 0), (1, 2), (2, 0))
+        e1 = b1.edge()
+        e2 = b2.edge()
+        assert not e1.geom_equal(e2)
+
+    def test_different_degree(self):
+        """Bezier curves with different degrees (different number of control points)."""
+        # Quadratic (degree 2, 3 points)
+        b1 = Bezier((0, 0), (1, 1), (2, 0))
+        # Cubic (degree 3, 4 points) - adjusted to have same endpoints
+        b2 = Bezier((0, 0), (0.5, 1), (1.5, 1), (2, 0))
+        e1 = b1.edge()
+        e2 = b2.edge()
+        assert e1.geom_type == GeomType.BEZIER
+        assert e2.geom_type == GeomType.BEZIER
+        assert not e1.geom_equal(e2)
+
+    def test_rational_bezier_different_weights(self):
+        """Rational Bezier curves with different weights."""
+        pts = [(0, 0, 0), (1, 1, 0), (2, 0, 0)]
+
+        # Create rational Bezier with weights [1, 2, 1]
+        e1 = Edge.make_bezier(*pts, weights=[1.0, 2.0, 1.0])
+
+        # Create rational Bezier with weights [1, 3, 1]
+        e2 = Edge.make_bezier(*pts, weights=[1.0, 3.0, 1.0])
+
+        assert e1.geom_type == GeomType.BEZIER
+        assert not e1.geom_equal(e2)
+
+
+class TestEdgeGeomEqualBSpline:
+    """Tests for Edge.geom_equal with BSPLINE type."""
+
+    def test_same_spline(self):
+        v = [Vertex(p) for p in ((-2, 0), (-1, 0), (0, 0), (1, 0), (2, 0))]
+        s1 = Spline(*v)
+        s2 = Spline(*v)
+        e1 = s1.edge()
+        e2 = s2.edge()
+        assert e1.geom_type == GeomType.BSPLINE
+        assert e1.geom_equal(e2)
+
+    def test_different_spline(self):
+        v1 = [Vertex(p) for p in ((-2, 0), (-1, 0), (0, 0), (1, 0), (2, 0))]
+        v2 = [Vertex(p) for p in ((-2, 0), (-1, 1), (0, 0), (1, 0), (2, 0))]
+        s1 = Spline(*v1)
+        s2 = Spline(*v2)
+        e1 = s1.edge()
+        e2 = s2.edge()
+        assert not e1.geom_equal(e2)
+
+    def test_complex_spline(self):
+        v = [
+            Vertex(p)
+            for p in (
+                (-2, 0),
+                (-1, 0),
+                (0, 0),
+                (1, 0),
+                (2, 0),
+                (3, 0.1),
+                (4, 1),
+                (5, 2.2),
+                (6, 3),
+                (7, 2),
+                (8, -1),
+            )
+        ]
+        s1 = Spline(*v)
+        s2 = Spline(*v)
+        e1 = s1.edge()
+        e2 = s2.edge()
+        assert e1.geom_equal(e2)
+
+    def test_different_periodicity(self):
+        """BSplines with different periodicity (periodic vs non-periodic)."""
+        # Same control points, different periodicity
+        pts = [(0, 0), (1, 1), (2, 0), (1, -1)]
+
+        e1 = Edge.make_spline(pts, periodic=False)
+        e2 = Edge.make_spline(pts, periodic=True)
+
+        assert e1.geom_type == GeomType.BSPLINE
+        assert e2.geom_type == GeomType.BSPLINE
+        # Different periodicity means not equal
+        assert not e1.geom_equal(e2)
+
+    def test_different_pole_count(self):
+        """BSplines with different number of poles."""
+        # 5 points
+        v1 = [Vertex(p) for p in ((0, 0), (1, 1), (2, 0), (3, 1), (4, 0))]
+        # 6 points with same endpoints
+        v2 = [
+            Vertex(p)
+            for p in ((0, 0), (0.8, 0.8), (1.6, 0.2), (2.4, 0.8), (3.2, 0.2), (4, 0))
+        ]
+        s1 = Spline(*v1)
+        s2 = Spline(*v2)
+        e1 = s1.edge()
+        e2 = s2.edge()
+        assert e1.geom_type == GeomType.BSPLINE
+        assert e2.geom_type == GeomType.BSPLINE
+        assert not e1.geom_equal(e2)
+
+    def test_different_knot_values(self):
+        """BSplines with different internal knot positions have different shapes."""
+        from OCP.Geom import Geom_BSplineCurve
+        from OCP.TColgp import TColgp_Array1OfPnt
+        from OCP.TColStd import TColStd_Array1OfReal, TColStd_Array1OfInteger
+        from OCP.gp import gp_Pnt
+        from OCP.BRepBuilderAPI import BRepBuilderAPI_MakeEdge
+
+        # 5 poles for degree 3 with one internal knot
+        poles = TColgp_Array1OfPnt(1, 5)
+        poles.SetValue(1, gp_Pnt(0, 0, 0))
+        poles.SetValue(2, gp_Pnt(1, 2, 0))
+        poles.SetValue(3, gp_Pnt(2, 2, 0))
+        poles.SetValue(4, gp_Pnt(3, 2, 0))
+        poles.SetValue(5, gp_Pnt(4, 0, 0))
+
+        mults = TColStd_Array1OfInteger(1, 3)
+        mults.SetValue(1, 4)
+        mults.SetValue(2, 1)  # Internal knot
+        mults.SetValue(3, 4)
+
+        # Internal knot at 0.5
+        knots1 = TColStd_Array1OfReal(1, 3)
+        knots1.SetValue(1, 0.0)
+        knots1.SetValue(2, 0.5)
+        knots1.SetValue(3, 1.0)
+        curve1 = Geom_BSplineCurve(poles, knots1, mults, 3, False)
+        e1 = Edge(BRepBuilderAPI_MakeEdge(curve1).Edge())
+
+        # Internal knot at 0.3 - different position changes shape!
+        knots2 = TColStd_Array1OfReal(1, 3)
+        knots2.SetValue(1, 0.0)
+        knots2.SetValue(2, 0.3)
+        knots2.SetValue(3, 1.0)
+        curve2 = Geom_BSplineCurve(poles, knots2, mults, 3, False)
+        e2 = Edge(BRepBuilderAPI_MakeEdge(curve2).Edge())
+
+        assert e1.geom_type == GeomType.BSPLINE
+        # Different internal knot position = different geometric shape
+        assert (e1 @ 0.5) != (e2 @ 0.5)
+        assert not e1.geom_equal(e2)
+
+    def test_different_multiplicities(self):
+        """BSplines with same poles/knots but different multiplicities have different shapes."""
+        from OCP.Geom import Geom_BSplineCurve
+        from OCP.TColgp import TColgp_Array1OfPnt
+        from OCP.TColStd import TColStd_Array1OfReal, TColStd_Array1OfInteger
+        from OCP.gp import gp_Pnt
+        from OCP.BRepBuilderAPI import BRepBuilderAPI_MakeEdge
+
+        # Same 7 poles for both curves
+        poles = TColgp_Array1OfPnt(1, 7)
+        poles.SetValue(1, gp_Pnt(0, 0, 0))
+        poles.SetValue(2, gp_Pnt(1, 2, 0))
+        poles.SetValue(3, gp_Pnt(2, 1, 0))
+        poles.SetValue(4, gp_Pnt(3, 2, 0))
+        poles.SetValue(5, gp_Pnt(4, 1, 0))
+        poles.SetValue(6, gp_Pnt(5, 2, 0))
+        poles.SetValue(7, gp_Pnt(6, 0, 0))
+
+        # Same 4 knots for both curves
+        knots = TColStd_Array1OfReal(1, 4)
+        knots.SetValue(1, 0.0)
+        knots.SetValue(2, 0.33)
+        knots.SetValue(3, 0.67)
+        knots.SetValue(4, 1.0)
+
+        # Multiplicities [4, 1, 2, 4] - sum = 11 = 7 + 3 + 1
+        mults1 = TColStd_Array1OfInteger(1, 4)
+        mults1.SetValue(1, 4)
+        mults1.SetValue(2, 1)
+        mults1.SetValue(3, 2)
+        mults1.SetValue(4, 4)
+        curve1 = Geom_BSplineCurve(poles, knots, mults1, 3, False)
+        e1 = Edge(BRepBuilderAPI_MakeEdge(curve1).Edge())
+
+        # Multiplicities [4, 2, 1, 4] - same sum, swapped internal mults
+        mults2 = TColStd_Array1OfInteger(1, 4)
+        mults2.SetValue(1, 4)
+        mults2.SetValue(2, 2)
+        mults2.SetValue(3, 1)
+        mults2.SetValue(4, 4)
+        curve2 = Geom_BSplineCurve(poles, knots, mults2, 3, False)
+        e2 = Edge(BRepBuilderAPI_MakeEdge(curve2).Edge())
+
+        assert e1.geom_type == GeomType.BSPLINE
+        assert e2.geom_type == GeomType.BSPLINE
+        # Same poles, same knots, different multiplicities = different shape
+        assert (e1 @ 0.5) != (e2 @ 0.5)
+        assert not e1.geom_equal(e2)
+
+    def test_rational_bspline_different_weights(self):
+        """Rational BSplines with different weights."""
+        from OCP.Geom import Geom_BSplineCurve
+        from OCP.TColgp import TColgp_Array1OfPnt
+        from OCP.TColStd import TColStd_Array1OfReal, TColStd_Array1OfInteger
+        from OCP.gp import gp_Pnt
+        from OCP.BRepBuilderAPI import BRepBuilderAPI_MakeEdge
+
+        poles = TColgp_Array1OfPnt(1, 4)
+        poles.SetValue(1, gp_Pnt(0, 0, 0))
+        poles.SetValue(2, gp_Pnt(1, 1, 0))
+        poles.SetValue(3, gp_Pnt(2, 1, 0))
+        poles.SetValue(4, gp_Pnt(3, 0, 0))
+
+        knots = TColStd_Array1OfReal(1, 2)
+        knots.SetValue(1, 0.0)
+        knots.SetValue(2, 1.0)
+        mults = TColStd_Array1OfInteger(1, 2)
+        mults.SetValue(1, 4)
+        mults.SetValue(2, 4)
+
+        # Weights [1, 2, 2, 1]
+        weights1 = TColStd_Array1OfReal(1, 4)
+        weights1.SetValue(1, 1.0)
+        weights1.SetValue(2, 2.0)
+        weights1.SetValue(3, 2.0)
+        weights1.SetValue(4, 1.0)
+        curve1 = Geom_BSplineCurve(poles, weights1, knots, mults, 3, False)
+        e1 = Edge(BRepBuilderAPI_MakeEdge(curve1).Edge())
+
+        # Weights [1, 3, 3, 1]
+        weights2 = TColStd_Array1OfReal(1, 4)
+        weights2.SetValue(1, 1.0)
+        weights2.SetValue(2, 3.0)
+        weights2.SetValue(3, 3.0)
+        weights2.SetValue(4, 1.0)
+        curve2 = Geom_BSplineCurve(poles, weights2, knots, mults, 3, False)
+        e2 = Edge(BRepBuilderAPI_MakeEdge(curve2).Edge())
+
+        assert e1.geom_type == GeomType.BSPLINE
+        assert not e1.geom_equal(e2)
+
+
+class TestEdgeGeomEqualOffset:
+    """Tests for Edge.geom_equal with OFFSET type."""
+
+    def test_same_offset(self):
+        v = [Vertex(p) for p in ((0, 0), (1, 1), (2, 0), (3, 1))]
+        s = Spline(*v)
+        w = Wire([s.edge()])
+        offset_wire1 = w.offset_2d(0.1)
+        offset_wire2 = w.offset_2d(0.1)
+
+        offset_edges1 = [
+            e for e in offset_wire1.edges() if e.geom_type == GeomType.OFFSET
+        ]
+        offset_edges2 = [
+            e for e in offset_wire2.edges() if e.geom_type == GeomType.OFFSET
+        ]
+
+        assert len(offset_edges1) > 0
+        assert offset_edges1[0].geom_equal(offset_edges2[0])
+
+    def test_different_offset_value(self):
+        v = [Vertex(p) for p in ((0, 0), (1, 1), (2, 0), (3, 1))]
+        s = Spline(*v)
+        w = Wire([s.edge()])
+        offset_wire1 = w.offset_2d(0.1)
+        offset_wire2 = w.offset_2d(0.2)
+
+        offset_edges1 = [
+            e for e in offset_wire1.edges() if e.geom_type == GeomType.OFFSET
+        ]
+        offset_edges2 = [
+            e for e in offset_wire2.edges() if e.geom_type == GeomType.OFFSET
+        ]
+
+        assert not offset_edges1[0].geom_equal(offset_edges2[0])
+
+    def test_different_offset_direction(self):
+        """Offset curves with different offset directions (on different planes)."""
+        from build123d import Axis
+
+        v = [Vertex(p) for p in ((0, 0), (1, 1), (2, 0), (3, 1))]
+        s = Spline(*v)
+        w = Wire([s.edge()])
+
+        # Offset on XY plane (Z direction)
+        offset_wire1 = w.offset_2d(0.1)
+        offset_edges1 = [
+            e for e in offset_wire1.edges() if e.geom_type == GeomType.OFFSET
+        ]
+
+        # Rotate wire 90 degrees around X axis to put it on XZ plane
+        w_rotated = w.rotate(Axis.X, 90)
+        offset_wire2 = w_rotated.offset_2d(0.1)
+        offset_edges2 = [
+            e for e in offset_wire2.edges() if e.geom_type == GeomType.OFFSET
+        ]
+
+        if len(offset_edges1) > 0 and len(offset_edges2) > 0:
+            # Different directions means not equal
+            assert not offset_edges1[0].geom_equal(offset_edges2[0])
+
+
+class TestEdgeGeomEqualTolerance:
+    """Tests for tolerance behavior in Edge.geom_equal."""
+
+    def test_circle_radius_within_tolerance(self):
+        """Circle radii differing by less than tolerance are equal."""
+        e1 = Edge.make_circle(10.0)
+        e2 = Edge.make_circle(10.0 + 1e-7)  # Within default tol=1e-6
+        assert e1.geom_equal(e2)
+
+    def test_circle_radius_outside_tolerance(self):
+        """Circle radii differing by more than tolerance are not equal."""
+        e1 = Edge.make_circle(10.0)
+        e2 = Edge.make_circle(10.0 + 1e-5)  # Outside default tol=1e-6
+        assert not e1.geom_equal(e2)
+
+    def test_tol_parameter_accepted(self):
+        """The tol parameter is accepted by geom_equal.
+
+        Note: The tol parameter affects property comparisons (radius, focal length,
+        weights, knots, offset values) but endpoint comparison always uses Vector's
+        built-in TOLERANCE (1e-6). Since most geometric differences also change
+        endpoints, the custom tol has limited practical effect.
+        """
+        e1 = Edge.make_line((0, 0), (1, 1))
+        e2 = Edge.make_line((0, 0), (1, 1))
+        # Parameter is accepted
+        assert e1.geom_equal(e2, tol=1e-9)
+        assert e1.geom_equal(e2, tol=0.1)
+
+    def test_line_endpoint_within_tolerance(self):
+        """Line endpoints differing by less than tolerance are equal."""
+        e1 = Edge.make_line((0, 0, 0), (1, 1, 1))
+        e2 = Edge.make_line((0, 0, 0), (1 + 1e-7, 1, 1))
+        assert e1.geom_equal(e2)
+
+    def test_line_endpoint_outside_tolerance(self):
+        """Line endpoints differing by more than tolerance are not equal."""
+        e1 = Edge.make_line((0, 0, 0), (1, 1, 1))
+        e2 = Edge.make_line((0, 0, 0), (1.001, 1, 1))
+        assert not e1.geom_equal(e2)
+
+
+class TestEdgeGeomEqualReversed:
+    """Tests for reversed edge comparison."""
+
+    def test_line_reversed_not_equal(self):
+        """Reversed line is not equal (different direction)."""
+        e1 = Edge.make_line((0, 0), (1, 1))
+        e2 = Edge.make_line((1, 1), (0, 0))
+        assert not e1.geom_equal(e2)
+
+    def test_arc_reversed_not_equal(self):
+        """Reversed arc is not equal."""
+        e1 = Edge.make_circle(10, start_angle=0, end_angle=90)
+        e2 = Edge.make_circle(10, start_angle=0, end_angle=90).reversed()
+        # Reversed edge has swapped start/end points
+        assert not e1.geom_equal(e2)
+
+    def test_spline_reversed_not_equal(self):
+        """Reversed spline is not equal."""
+        pts = [(0, 0), (1, 1), (2, 0), (3, 1)]
+        s = Spline(*[Vertex(p) for p in pts])
+        e1 = s.edge()
+        e2 = e1.reversed()
+        assert not e1.geom_equal(e2)
+
+
+class TestEdgeGeomEqualArcVariations:
+    """Tests for arc edge cases."""
+
+    def test_full_circle_equal(self):
+        """Two full circles are equal."""
+        e1 = Edge.make_circle(10)
+        e2 = Edge.make_circle(10)
+        assert e1.geom_equal(e2)
+
+    def test_arc_different_start_same_sweep(self):
+        """Arcs with different start angles but same sweep are not equal."""
+        e1 = Edge.make_circle(10, start_angle=0, end_angle=90)
+        e2 = Edge.make_circle(10, start_angle=90, end_angle=180)
+        # Same radius and sweep angle, but different positions
+        assert not e1.geom_equal(e2)
+
+    def test_arc_same_endpoints_different_direction(self):
+        """Arcs with same endpoints but opposite sweep direction."""
+        from build123d import AngularDirection
+
+        e1 = Edge.make_circle(
+            10,
+            start_angle=0,
+            end_angle=90,
+            angular_direction=AngularDirection.COUNTER_CLOCKWISE,
+        )
+        e2 = Edge.make_circle(
+            10,
+            start_angle=90,
+            end_angle=0,
+            angular_direction=AngularDirection.CLOCKWISE,
+        )
+        # These trace different paths (short arc vs long arc)
+        assert not e1.geom_equal(e2)
+
+
+class TestEdgeGeomEqualNumerical:
+    """Tests for numerical edge cases."""
+
+    def test_very_small_edge(self):
+        """Very small edges can be compared."""
+        e1 = Edge.make_line((0, 0), (1e-6, 1e-6))
+        e2 = Edge.make_line((0, 0), (1e-6, 1e-6))
+        assert e1.geom_equal(e2)
+
+    def test_very_large_coordinates(self):
+        """Edges with large coordinates can be compared."""
+        e1 = Edge.make_line((1e6, 1e6), (1e6 + 1, 1e6 + 1))
+        e2 = Edge.make_line((1e6, 1e6), (1e6 + 1, 1e6 + 1))
+        assert e1.geom_equal(e2)
+
+    def test_large_coordinates_small_difference(self):
+        """Small differences at large coordinates."""
+        e1 = Edge.make_line((1e6, 1e6), (1e6 + 1, 1e6 + 1))
+        e2 = Edge.make_line((1e6, 1e6), (1e6 + 1 + 1e-5, 1e6 + 1))
+        # Difference is above tolerance
+        assert not e1.geom_equal(e2)
+
+
+class TestEdgeGeomEqual3DPositioning:
+    """Tests for 3D positioning edge cases."""
+
+    def test_same_shape_different_z(self):
+        """Same 2D shape at different Z levels are not equal."""
+        e1 = Edge.make_circle(10)
+        e2 = Edge.make_circle(10).locate(Location((0, 0, 5)))
+        assert not e1.geom_equal(e2)
+
+    def test_line_in_3d(self):
+        """3D lines with same geometry are equal."""
+        e1 = Edge.make_line((0, 0, 0), (1, 2, 3))
+        e2 = Edge.make_line((0, 0, 0), (1, 2, 3))
+        assert e1.geom_equal(e2)
+
+    def test_line_in_3d_different(self):
+        """3D lines with different Z are not equal."""
+        e1 = Edge.make_line((0, 0, 0), (1, 1, 0))
+        e2 = Edge.make_line((0, 0, 0), (1, 1, 1))
+        assert not e1.geom_equal(e2)
+
+
+class TestEdgeGeomEqualSplineVariations:
+    """Tests for BSpline edge cases."""
+
+    def test_spline_control_point_within_tolerance(self):
+        """Splines with control points within tolerance are equal."""
+        pts1 = [(0, 0), (1, 1), (2, 0)]
+        pts2 = [(0, 0), (1 + 1e-7, 1), (2, 0)]
+        e1 = Spline(*[Vertex(p) for p in pts1]).edge()
+        e2 = Spline(*[Vertex(p) for p in pts2]).edge()
+        assert e1.geom_equal(e2)
+
+    def test_spline_control_point_outside_tolerance(self):
+        """Splines with control points outside tolerance are not equal."""
+        pts1 = [(0, 0), (1, 1), (2, 0)]
+        pts2 = [(0, 0), (1.001, 1), (2, 0)]
+        e1 = Spline(*[Vertex(p) for p in pts1]).edge()
+        e2 = Spline(*[Vertex(p) for p in pts2]).edge()
+        assert not e1.geom_equal(e2)
+
+    def test_spline_different_point_count(self):
+        """Splines with different number of control points are not equal."""
+        pts1 = [(0, 0), (1, 1), (2, 0)]
+        pts2 = [(0, 0), (0.5, 0.5), (1, 1), (2, 0)]
+        e1 = Spline(*[Vertex(p) for p in pts1]).edge()
+        e2 = Spline(*[Vertex(p) for p in pts2]).edge()
+        # Different number of poles
+        assert not e1.geom_equal(e2)
+
+
+class TestEdgeGeomEqualUnknownType:
+    """Tests for the fallback case (OTHER/unknown geom types)."""
+
+    def test_interpolation_points_used(self):
+        """For unknown types, sample points are compared."""
+        # Create edges that would use the fallback path
+        # Most common types are handled, but we can test the parameter
+        e1 = Edge.make_line((0, 0), (1, 1))
+        e2 = Edge.make_line((0, 0), (1, 1))
+        # Even with different num_interpolation_points, these should be equal
+        assert e1.geom_equal(e2, num_interpolation_points=3)
+        assert e1.geom_equal(e2, num_interpolation_points=10)
+
+
+class TestWireGeomEqual:
+    """Tests for Wire.geom_equal method."""
+
+    def test_same_rectangle_wire(self):
+        r1 = Rectangle(10, 5)
+        r2 = Rectangle(10, 5)
+        assert r1.wire().geom_equal(r2.wire())
+
+    def test_different_rectangle_wire(self):
+        r1 = Rectangle(10, 5)
+        r2 = Rectangle(10, 6)
+        assert not r1.wire().geom_equal(r2.wire())
+
+    def test_same_spline_wire(self):
+        v = [Vertex(p) for p in ((0, 0), (1, 1), (2, 0), (3, 1))]
+        s1 = Spline(*v)
+        s2 = Spline(*v)
+        w1 = Wire([s1.edge()])
+        w2 = Wire([s2.edge()])
+        assert w1.geom_equal(w2)
+
+    def test_different_edge_count(self):
+        r1 = Rectangle(10, 5)
+        e = Edge.make_line((0, 0), (10, 0))
+        w1 = r1.wire()
+        w2 = Wire([e])
+        assert not w1.geom_equal(w2)
+
+    def test_identical_edge_objects(self):
+        """Two wires sharing the same edge objects."""
+        e1 = Edge.make_line((0, 0), (1, 0))
+        e2 = Edge.make_line((1, 0), (1, 1))
+        e3 = Edge.make_line((1, 1), (0, 0))
+        w1 = Wire([e1, e2, e3])
+        w2 = Wire([e1, e2, e3])  # Same edge objects
+        assert w1.geom_equal(w2)
+
+    def test_geometrically_equal_edges(self):
+        """Two wires with geometrically equal but distinct edge objects."""
+        # Wire 1
+        e1a = Edge.make_line((0, 0), (1, 0))
+        e2a = Edge.make_line((1, 0), (1, 1))
+        e3a = Edge.make_line((1, 1), (0, 0))
+        w1 = Wire([e1a, e2a, e3a])
+        # Wire 2 - same geometry, different objects
+        e1b = Edge.make_line((0, 0), (1, 0))
+        e2b = Edge.make_line((1, 0), (1, 1))
+        e3b = Edge.make_line((1, 1), (0, 0))
+        w2 = Wire([e1b, e2b, e3b])
+        assert w1.geom_equal(w2)
+
+    def test_edges_different_start_point(self):
+        """Two closed wires with same geometry but different starting vertex are not equal."""
+        # Wire 1: starts at (0,0)
+        e1a = Edge.make_line((0, 0), (1, 0))
+        e2a = Edge.make_line((1, 0), (1, 1))
+        e3a = Edge.make_line((1, 1), (0, 0))
+        w1 = Wire([e1a, e2a, e3a])
+        # Wire 2: starts at (1,1) due to different edge order in constructor
+        e3b = Edge.make_line((1, 1), (0, 0))
+        e1b = Edge.make_line((0, 0), (1, 0))
+        e2b = Edge.make_line((1, 0), (1, 1))
+        w2 = Wire([e3b, e1b, e2b])
+        # Different starting point means not equal
+        assert not w1.geom_equal(w2)
+
+    def test_one_edge_reversed(self):
+        """Two wires where one has an edge with reversed direction."""
+        # Wire 1: all edges in forward direction
+        e1a = Edge.make_line((0, 0), (1, 0))
+        e2a = Edge.make_line((1, 0), (1, 1))
+        e3a = Edge.make_line((1, 1), (0, 0))
+        w1 = Wire([e1a, e2a, e3a])
+        # Wire 2: middle edge is reversed (direction (1,1) -> (1,0) instead of (1,0) -> (1,1))
+        e1b = Edge.make_line((0, 0), (1, 0))
+        e2b = Edge.make_line((1, 1), (1, 0))  # Reversed!
+        e3b = Edge.make_line((1, 1), (0, 0))
+        w2 = Wire([e1b, e2b, e3b])
+        # order_edges should correct the orientation
+        assert w1.geom_equal(w2)
+
+    def test_closed_wire(self):
+        """Two closed wires with same geometry."""
+        w1 = Wire(
+            [
+                Edge.make_line((0, 0), (2, 0)),
+                Edge.make_line((2, 0), (2, 2)),
+                Edge.make_line((2, 2), (0, 2)),
+                Edge.make_line((0, 2), (0, 0)),
+            ]
+        )
+        w2 = Wire(
+            [
+                Edge.make_line((0, 0), (2, 0)),
+                Edge.make_line((2, 0), (2, 2)),
+                Edge.make_line((2, 2), (0, 2)),
+                Edge.make_line((0, 2), (0, 0)),
+            ]
+        )
+        assert w1.is_closed
+        assert w2.is_closed
+        assert w1.geom_equal(w2)
+
+    def test_mixed_edge_types(self):
+        """Wires with mixed edge types (lines and arcs)."""
+        # Wire with line + arc + line
+        e1a = Edge.make_line((0, 0), (1, 0))
+        e2a = Edge.make_circle(0.5, start_angle=0, end_angle=180).locate(
+            Location((1.5, 0, 0))
+        )
+        e3a = Edge.make_line((2, 0), (3, 0))
+        w1 = Wire([e1a, e2a, e3a])
+
+        e1b = Edge.make_line((0, 0), (1, 0))
+        e2b = Edge.make_circle(0.5, start_angle=0, end_angle=180).locate(
+            Location((1.5, 0, 0))
+        )
+        e3b = Edge.make_line((2, 0), (3, 0))
+        w2 = Wire([e1b, e2b, e3b])
+
+        assert w1.geom_equal(w2)
+
+    def test_mixed_edge_types_different(self):
+        """Wires with mixed edge types that differ."""
+        # Wire 1: line + arc
+        e1a = Edge.make_line((0, 0), (1, 0))
+        e2a = Edge.make_circle(0.5, start_angle=0, end_angle=180).locate(
+            Location((1.5, 0, 0))
+        )
+        w1 = Wire([e1a, e2a])
+
+        # Wire 2: line + different arc (different radius)
+        e1b = Edge.make_line((0, 0), (1, 0))
+        e2b = Edge.make_circle(0.6, start_angle=0, end_angle=180).locate(
+            Location((1.6, 0, 0))
+        )
+        w2 = Wire([e1b, e2b])
+
+        assert not w1.geom_equal(w2)
+
+    def test_all_edges_reversed_not_equal(self):
+        """Wire traced in opposite direction is not equal."""
+        # Wire 1: (0,0) -> (3,0)
+        e1a = Edge.make_line((0, 0), (1, 0))
+        e2a = Edge.make_line((1, 0), (2, 1))
+        e3a = Edge.make_line((2, 1), (3, 0))
+        w1 = Wire([e1a, e2a, e3a])
+
+        # Wire 2: (3,0) -> (0,0) - same path but opposite direction
+        e1b = Edge.make_line((1, 0), (0, 0))
+        e2b = Edge.make_line((2, 1), (1, 0))
+        e3b = Edge.make_line((3, 0), (2, 1))
+        w2 = Wire([e3b, e2b, e1b])
+
+        assert not w1.geom_equal(w2)
+
+    def test_open_wire_different_start(self):
+        """Open wires with same edges but different starting edge - should not match."""
+        # For open wires, the start matters
+        e1 = Edge.make_line((0, 0), (1, 0))
+        e2 = Edge.make_line((1, 0), (2, 1))
+        e3 = Edge.make_line((2, 1), (3, 0))
+        w1 = Wire([e1, e2, e3])
+
+        # Different edges entirely (shifted)
+        e4 = Edge.make_line((1, 0), (2, 0))
+        e5 = Edge.make_line((2, 0), (3, 1))
+        e6 = Edge.make_line((3, 1), (4, 0))
+        w2 = Wire([e4, e5, e6])
+
+        assert not w1.geom_equal(w2)
+
+    def test_wire_with_spline_edges(self):
+        """Wires containing spline edges."""
+        pts1 = [(0, 0), (1, 1), (2, 0), (3, 1), (4, 0)]
+        pts2 = [(4, 0), (5, 1), (6, 0)]
+
+        s1a = Spline(*[Vertex(p) for p in pts1])
+        s2a = Spline(*[Vertex(p) for p in pts2])
+        w1 = Wire([s1a.edge(), s2a.edge()])
+
+        s1b = Spline(*[Vertex(p) for p in pts1])
+        s2b = Spline(*[Vertex(p) for p in pts2])
+        w2 = Wire([s1b.edge(), s2b.edge()])
+
+        assert w1.geom_equal(w2)
+
+    def test_single_edge_wire(self):
+        """Wires with single edge."""
+        w1 = Wire([Edge.make_line((0, 0), (5, 5))])
+        w2 = Wire([Edge.make_line((0, 0), (5, 5))])
+        assert w1.geom_equal(w2)
+
+    def test_single_edge_wire_reversed_not_equal(self):
+        """Single edge wire vs reversed single edge wire are not equal."""
+        w1 = Wire([Edge.make_line((0, 0), (5, 5))])
+        w2 = Wire([Edge.make_line((5, 5), (0, 0))])
+        # Opposite direction means not equal
+        assert not w1.geom_equal(w2)
+
+
+class TestGeomEqualTypeMismatch:
+    """Tests for type mismatch cases."""
+
+    def test_edge_vs_non_edge(self):
+        e = Edge.make_line((0, 0), (1, 1))
+        w = Wire([e])
+        # Edge.geom_equal should return False for non-Edge
+        assert not e.geom_equal(w)
+
+    def test_wire_vs_non_wire(self):
+        e = Edge.make_line((0, 0), (1, 1))
+        w = Wire([e])
+        # Wire.geom_equal should return False for non-Wire
+        assert not w.geom_equal(e)
+
+    def test_different_geom_types(self):
+        line = Edge.make_line((0, 0, 0), (1, 1, 1))
+        circle = Circle(10).edge()
+        assert not line.geom_equal(circle)

tests/test_direct_api/test_intersection.py

@@ -16,14 +16,19 @@ class Case:
     expected: list | Vector | Location | Axis | Plane
     name: str
     xfail: None | str = None
+    include_touched: bool = False
 
 
 @pytest.mark.skip
-def run_test(obj, target, expected):
+def run_test(obj, target, expected, include_touched=False):
+    # Only Shape objects support include_touched parameter
+    kwargs = {}
+    if include_touched and isinstance(obj, Shape):
+        kwargs["include_touched"] = include_touched
     if isinstance(target, list):
-        result = obj.intersect(*target)
+        result = obj.intersect(*target, **kwargs)
     else:
-        result = obj.intersect(target)
+        result = obj.intersect(target, **kwargs)
     if INTERSECT_DEBUG:
         show([obj, target, result])
     if expected is None:
@@ -50,11 +55,15 @@ def make_params(matrix):
             marks = [pytest.mark.xfail(reason=case.xfail)]
         else:
             marks = []
-        uid = f"{i} {obj_type}, {tar_type}, {case.name}"
-        params.append(pytest.param(case.object, case.target, case.expected, marks=marks, id=uid))
-        if tar_type != obj_type and not isinstance(case.target, list):
-            uid = f"{i + 1} {tar_type}, {obj_type}, {case.name}"
-            params.append(pytest.param(case.target, case.object, case.expected, marks=marks, id=uid))
+        # Add include_touched info to test id if specified
+        touched_suffix = ", touched" if case.include_touched else ""
+        uid = f"{i} {obj_type}, {tar_type}, {case.name}{touched_suffix}"
+        params.append(pytest.param(case.object, case.target, case.expected, case.include_touched, marks=marks, id=uid))
+        # Swap obj and target to test symmetry, but NOT for include_touched tests
+        # (swapping may change behavior with boundary contacts)
+        if tar_type != obj_type and not isinstance(case.target, list) and not case.include_touched:
+            uid = f"{i + 1} {tar_type}, {obj_type}, {case.name}{touched_suffix}"
+            params.append(pytest.param(case.target, case.object, case.expected, case.include_touched, marks=marks, id=uid))
 
     return params
 
@@ -118,9 +127,9 @@ geometry_matrix = [
     Case(lc1, lc1, Location, "coincident, co-z", None),
 ]
 
-@pytest.mark.parametrize("obj, target, expected", make_params(geometry_matrix))
-def test_geometry(obj, target, expected):
-    run_test(obj, target, expected)
+@pytest.mark.parametrize("obj, target, expected, include_touched", make_params(geometry_matrix))
+def test_geometry(obj, target, expected, include_touched):
+    run_test(obj, target, expected, include_touched)
 
 
 # Shape test matrices
@@ -147,9 +156,9 @@ shape_0d_matrix = [
     Case(vt1, [vt1, lc1], [Vertex], "multi to_intersect, coincident", None),
 ]
 
-@pytest.mark.parametrize("obj, target, expected", make_params(shape_0d_matrix))
-def test_shape_0d(obj, target, expected):
-    run_test(obj, target, expected)
+@pytest.mark.parametrize("obj, target, expected, include_touched", make_params(shape_0d_matrix))
+def test_shape_0d(obj, target, expected, include_touched):
+    run_test(obj, target, expected, include_touched)
 
 
 # 1d Shapes
@@ -216,9 +225,9 @@ shape_1d_matrix = [
     Case(wi5, [ed1, Vector(1, 0)], [Vertex], "multi to_intersect, multi intersect", None),
 ]
 
-@pytest.mark.parametrize("obj, target, expected", make_params(shape_1d_matrix))
-def test_shape_1d(obj, target, expected):
-    run_test(obj, target, expected)
+@pytest.mark.parametrize("obj, target, expected, include_touched", make_params(shape_1d_matrix))
+def test_shape_1d(obj, target, expected, include_touched):
+    run_test(obj, target, expected, include_touched)
 
 
 # 2d Shapes
@@ -229,6 +238,16 @@ fc4 = Rot(Z=45) * Rectangle(5, 5).face()
 fc5 = Pos(2.5, 2.5, 2.5) * Rot(0, 90) * Rectangle(5, 5).face()
 fc6 = Pos(2.5, 2.5) * Rot(0, 90, 45, Extrinsic.XYZ) * Rectangle(5, 5).face()
 fc7 = (Rot(90) * Cylinder(2, 4)).faces().filter_by(GeomType.CYLINDER)[0]
+fc8 = make_face(
+    Polyline(
+        (-1.5, 1, 1),
+        (-1.5, -1, 1),
+        (3.5, -1, -1),
+        (3.5, 1, -1),
+        (-1.5, 1, 1),
+    )
+)
+fc9 = Pos(-2) * mirror(fc8, Plane.XY)
 
 fc11 = Rectangle(4, 4).face()
 fc22 = sweep(Rot(90) * CenterArc((0, 0), 2, 0, 180), Line((0, 2), (0, -2)))
@@ -237,6 +256,13 @@ sh2 = Pos(Z=1) * sh1
 sh3 = Shell([Pos(-4) * fc11, fc22, Pos(2, 0, -2) * Rot(0, 90) * fc11])
 sh4 = Shell([Pos(-4) * fc11, fc22, Pos(4) * fc11])
 sh5 = Pos(Z=1) * Shell([Pos(-2, 0, -2) * Rot(0, -90) * fc11, fc22, Pos(2, 0, -2) * Rot(0, 90) * fc11])
+sh6 = Box(2, 2, 2).shell()
+
+# Shell tangent touch test objects (half spheres)
+_half_sphere_solid = Sphere(1) & Pos(0, 0, 0.5) * Box(3, 3, 2)
+sh7 = Shell(_half_sphere_solid.faces())
+sh8 = Pos(2, 0, 0) * sh7  # tangent at (1, 0, 0)
+fc10 = Pos(1, 0, 0) * (Rot(0, 90, 0) * Rectangle(2, 2).face())  # tangent to sphere at x=1
 
 shape_2d_matrix = [
     Case(fc1, vl2, None, "non-coincident", None),
@@ -272,7 +298,9 @@ shape_2d_matrix = [
     Case(fc1, fc3, [Edge], "intersecting", None),
     Case(fc1, fc4, [Face], "coplanar", None),
     Case(fc1, fc5, [Edge], "intersecting edge", None),
-    Case(fc1, fc6, [Vertex], "intersecting vertex", None),
+    # Face + Face crossing vertex: now requires include_touched
+    Case(fc1, fc6, None, "crossing vertex", None),
+    Case(fc1, fc6, [Vertex], "crossing vertex", None, True),
     Case(fc1, fc7, [Edge, Edge], "multi-intersecting", None),
     Case(fc7, Pos(Y=2) * fc7, [Face], "cyl intersecting", None),
 
@@ -281,20 +309,38 @@ shape_2d_matrix = [
     Case(sh1, fc1, [Face, Edge], "coplanar + intersecting", None),
     Case(sh4, fc1, [Face, Face], "2 coplanar", None),
     Case(sh5, fc1, [Edge, Edge], "2 intersecting", None),
+    Case(sh6, Pos(0,0,1) * fc1, [Face], "2 intersecting boundary", None),
+    Case(sh6, Pos(2, 1, 1) * sh6, [Face], "2 intersecting boundary", None),
+
+    # Shell + Face tangent touch
+    Case(sh7, fc10, None, "tangent touch", None),
+    Case(sh7, fc10, [Vertex], "tangent touch", None, True),
+    # Shell + Shell tangent touch
+    Case(sh7, sh8, None, "tangent touch", None),
+    Case(sh7, sh8, [Vertex], "tangent touch", None, True),
 
     Case(fc1, [fc4, Pos(2, 2) * fc1], [Face], "multi to_intersect, intersecting", None),
     Case(fc1, [ed1, Pos(2.5, 2.5) * fc1], [Edge], "multi to_intersect, intersecting", None),
     Case(fc7, [wi5, fc1], [Vertex], "multi to_intersect, intersecting", None),
 ]
 
-@pytest.mark.parametrize("obj, target, expected", make_params(shape_2d_matrix))
-def test_shape_2d(obj, target, expected):
-    run_test(obj, target, expected)
+@pytest.mark.parametrize("obj, target, expected, include_touched", make_params(shape_2d_matrix))
+def test_shape_2d(obj, target, expected, include_touched):
+    run_test(obj, target, expected, include_touched)
 
 # 3d Shapes
 sl1 = Box(2, 2, 2).solid()
 sl2 = Pos(Z=5) * Box(2, 2, 2).solid()
 sl3 = Cylinder(2, 1).solid() - Cylinder(1.5, 1).solid()
+sl4 = Box(3, 1, 1)
+# T-shaped solid (box + thin plate) for testing coplanar face touches
+sl5 = Pos(0.5, 0, 1) * Box(1, 1, 1) + Pos(0.5, 0, 1) * Box(2, 0.1, 1)
+sl6 = Pos(2, 0, 1.5) * Box(2, 2, 1)
+# Overlapping boxes where coplanar face is part of intersection (not touch)
+sl7 = Pos(0, 0.1, 0) * Box(2, 2, 2)
+sl8 = Pos(1, 0, -1) * Box(4, 2, 1)
+# Extended T-shaped solid for testing coplanar edge touches
+sl9 = Box(2, 2, 2) + sl5
 
 wi7 = Wire([l1 := sl3.faces().sort_by(Axis.Z)[-1].edges()[0].trim(.3, .4),
           l2 := l1.trim(2, 3),
@@ -316,6 +362,7 @@ shape_3d_matrix = [
 
     Case(sl1, pl3, None, "non-coincident", None),
     Case(sl1, pl2, [Face], "intersecting", None),
+    Case(sl1, pl2.offset(1), [Face], "intersecting boondary", None),
 
     Case(sl2, vt1, None, "non-coincident", None),
     Case(Pos(2) * sl1, vt1, [Vertex], "contained", None),
@@ -323,42 +370,85 @@ shape_3d_matrix = [
 
     Case(sl1, ed3, None, "non-coincident", None),
     Case(sl1, ed1, [Edge], "intersecting", None),
-    Case(sl1, Pos(0, 1, 1) * ed1, [Edge], "edge collinear", "duplicate edges, BRepAlgoAPI_Common and _Section both return edge"),
-    Case(sl1, Pos(1, 1, 1) * ed1, [Vertex], "corner coincident", None),
+    Case(sl1, Pos(0, 1, 1) * ed1, [Edge], "edge collinear", None),  # xfail removed
+    # Solid + Edge corner coincident: now requires include_touched
+    Case(sl1, Pos(1, 1, 1) * ed1, None, "corner coincident", None),
+    Case(sl1, Pos(1, 1, 1) * ed1, [Vertex], "corner coincident", None, True),
     Case(Pos(2.1, 1) * sl1, ed4, [Edge, Edge], "multi-intersect", None),
+    Case(Pos(2.1, 1, -1) * sl1, ed4, [Edge, Edge], "multi-intersect, boundary", None),
 
     Case(Pos(2, .5, -1) * sl1, wi6, None, "non-coincident", None),
     Case(Pos(2, .5, 1) * sl1, wi6, [Edge, Edge], "multi-intersecting", None),
+    Case(Pos(2, .5, 2) * sl1, wi6, [Edge, Edge], "multi-intersecting, boundary", None),
     Case(sl3, wi7, [Edge, Edge], "multi-coincident, is_equal check", None),
 
     Case(sl2, fc1, None, "non-coincident", None),
     Case(sl1, fc1, [Face], "intersecting", None),
-    Case(Pos(3.5, 0, 1) * sl1, fc1, [Edge], "edge collinear", None),
-    Case(Pos(3.5, 3.5) * sl1, fc1, [Vertex], "corner coincident", None),
+    Case(Pos(0,0,-1) * sl1, fc1, [Face], "intersecting, boundary", None),
+    Case(Pos(0,0,1) * sl1, fc1, [Face], "intersecting, boundary", None),
+    # Solid + Face edge collinear: now requires include_touched
+    Case(Pos(3.5, 0, 1) * sl1, fc1, None, "edge collinear", None),
+    Case(Pos(3.5, 0, 1) * sl1, fc1, [Edge], "edge collinear", None, True),
+    # Solid + Face corner coincident: now requires include_touched
+    Case(Pos(3.5, 3.5) * sl1, fc1, None, "corner coincident", None),
+    Case(Pos(3.5, 3.5) * sl1, fc1, [Vertex], "corner coincident", None, True),
     Case(Pos(.9) * sl1, fc7, [Face, Face], "multi-intersecting", None),
+    Case(Pos(.9,1) * sl1, fc7, [Face, Face], "multi-intersecting", None),
+    Case(Pos(.9,1.5) * sl1, fc7, [Face, Face], "multi-intersecting", None),
 
     Case(sl2, sh1, None, "non-coincident", None),
     Case(Pos(-2) * sl1, sh1, [Face, Face], "multi-intersecting", None),
+    Case(Pos(-2) * sl1, sh1, [Face, Face], "multi-intersecting", None),
+    Case(Pos(-2,3) * sl1, sh1, None, "multi-intersecting", None),
+    Case(Pos(-2,3) * sl1, sh1, [Edge, Edge], "multi-intersecting", None, True),
 
     Case(sl1, sl2, None, "non-coincident", None),
     Case(sl1, Pos(1, 1, 1) * sl1, [Solid], "intersecting", None),
-    Case(sl1, Pos(2, 2, 1) * sl1, [Edge], "edge collinear", None),
-    Case(sl1, Pos(2, 2, 2) * sl1, [Vertex], "corner coincident", None),
+    # Solid + Solid edge collinear: now requires include_touched
+    Case(sl1, Pos(2, 2, 1) * sl1, None, "edge collinear", None),
+    Case(sl1, Pos(2, 2, 1) * sl1, [Edge], "edge collinear", None, True),
+    # Solid + Solid face collinear: now requires include_touched
+    Case(sl1, Pos(2, 1.5, 1) * sl1, None, "edge collinear", None),
+    Case(sl1, Pos(2, 1.5, 1) * sl1, [Face], "edge collinear", None, True),
+    # Solid + Solid corner coincident: now requires include_touched
+    Case(sl1, Pos(2, 2, 2) * sl1, None, "corner coincident", None),
+    Case(sl1, Pos(2, 2, 2) * sl1, [Vertex], "corner coincident", None, True),
     Case(sl1, Pos(.45) * sl3, [Solid, Solid], "multi-intersect", None),
+    # New test: Solid + Solid face coincident (touch)
+    Case(sl1, Pos(2, 0, 0) * sl1, None, "face coincident", None),
+    Case(sl1, Pos(2, 0, 0) * sl1, [Face], "face coincident", None, True),
 
     Case(Pos(1.5, 1.5) * sl1, [sl3, Pos(.5, .5) * sl1], [Solid], "multi to_intersect, intersecting", None),
+    Case(Pos(0, 1.5) * sl1, [sl3, Pos(.5, .5) * sl1], [Solid, Solid], "multi to_intersect, intersecting", None),
+    Case(Pos(0.5, 1.5) * sl1, [sl3, Pos(.5, .5) * sl1], [Solid, Solid], "multi to_intersect, intersecting", None),
+    Case(Pos(0.5, 1.5) * sl1, [sl3, Pos(.5, .5) * sl1], [Solid, Solid], "multi to_intersect, intersecting", None, True),
+
     Case(Pos(1.5, 1.5) * sl1, [sl3, Pos(Z=.5) * fc1], [Face], "multi to_intersect, intersecting", None),
+
+    # T-shaped solid with coplanar face touches (edges should be filtered)
+    Case(sl5, sl6, [Solid], "coplanar face touch", None),
+    Case(sl5, sl6, [Solid, Face, Face], "coplanar face touch", None, True),
+
+    # Overlapping boxes: coplanar face is part of intersection, not touch
+    Case(sl7, sl8, [Solid], "coplanar face filtered", None),
+    Case(sl7, sl8, [Solid], "coplanar face filtered", None, True),
+
+    # Extended T-shaped solid with coplanar edge touches
+    Case(sl9, sl6, [Solid], "coplanar edge touch", None),
+    Case(sl9, sl6, [Solid, Face, Face, Edge, Edge], "coplanar edge touch", None, True),
 ]
 
-@pytest.mark.parametrize("obj, target, expected", make_params(shape_3d_matrix))
-def test_shape_3d(obj, target, expected):
-    run_test(obj, target, expected)
+@pytest.mark.parametrize("obj, target, expected, include_touched", make_params(shape_3d_matrix))
+def test_shape_3d(obj, target, expected, include_touched):
+    run_test(obj, target, expected, include_touched)
 
 # Compound Shapes
 cp1 = Compound(GridLocations(5, 0, 2, 1) * Vertex())
 cp2 = Compound(GridLocations(5, 0, 2, 1) * Line((0, -1), (0, 1)))
 cp3 = Compound(GridLocations(5, 0, 2, 1) * Rectangle(2, 2))
 cp4 = Compound(GridLocations(5, 0, 2, 1) * Box(2, 2, 2))
+cp5 = Compound([fc8, fc9])
+cp6 = Compound(GridLocations(4, 0, 2, 1) * Rectangle(2, 2))
 
 cv1 = Curve() + [ed1, ed2, ed3]
 sk1 = Sketch() + [fc1, fc2, fc3]
@@ -367,48 +457,61 @@ pt1 = Part() + [sl1, sl2, sl3]
 
 shape_compound_matrix = [
     Case(cp1, vl1, None, "non-coincident", None),
-    Case(Pos(-.5) * cp1, vl1, [Vertex], "intersecting", None),
-
+    Case(Pos(-0.5) * cp1, vl1, [Vertex], "intersecting", None),
     Case(cp2, lc1, None, "non-coincident", None),
-    Case(Pos(-.5) * cp2, lc1, [Vertex], "intersecting", None),
-
+    Case(Pos(-0.5) * cp2, lc1, [Vertex], "intersecting", None),
     Case(Pos(Z=1) * cp3, ax1, None, "non-coincident", None),
     Case(cp3, ax1, [Edge, Edge], "intersecting", None),
-
     Case(Pos(Z=3) * cp4, pl2, None, "non-coincident", None),
     Case(cp4, pl2, [Face, Face], "intersecting", None),
-
+    Case(Pos(Z=1) * cp4, pl2, [Face, Face], "non-coincident, boundary", None),
+    Case(Pos(Z=-1) * cp4, pl2, [Face, Face], "non-coincident, boundary", None),
     Case(cp1, vt1, None, "non-coincident", None),
-    Case(Pos(-.5) * cp1, vt1, [Vertex], "intersecting", None),
-
+    Case(Pos(-0.5) * cp1, vt1, [Vertex], "intersecting", None),
     Case(Pos(Z=1) * cp2, ed1, None, "non-coincident", None),
     Case(cp2, ed1, [Vertex], "intersecting", None),
-
     Case(Pos(Z=1) * cp3, fc1, None, "non-coincident", None),
     Case(cp3, fc1, [Face, Face], "intersecting", None),
-
+    Case(Pos(1) * cp3, fc1, [Face, Edge], "intersectingPos(0.5), ", None),
     Case(Pos(Z=5) * cp4, sl1, None, "non-coincident", None),
     Case(Pos(2) * cp4, sl1, [Solid], "intersecting", None),
-
+    Case(cp4, sl4, None, "intersecting", None),
+    Case(cp4, sl4, [Face, Face], "intersecting", None, True),
+    Case(cp4, Pos(0, 1, 1) * sl4, [Face, Face], "intersecting", None, True),
+    Case(cp4, Pos(0, 1, 1.5) * sl4, [Edge, Edge], "intersecting", None, True),
+    Case(cp4, Pos(0, 1.5, 1.5) * sl4, [Vertex, Vertex], "intersecting", None, True),
     Case(cp1, Pos(Z=1) * cp1, None, "non-coincident", None),
     Case(cp1, cp2, [Vertex, Vertex], "intersecting", None),
     Case(cp2, cp3, [Edge, Edge], "intersecting", None),
+    Case(Pos(0, 2, 0) * cp2, cp3, [Vertex, Vertex], "intersecting", None),
     Case(cp3, cp4, [Face, Face], "intersecting", None),
-
-    Case(cp1, Compound(children=cp1.get_type(Vertex)), [Vertex, Vertex], "mixed child type", None),
-    Case(cp4, Compound(children=cp3.get_type(Face)), [Face, Face], "mixed child type", None),
-
+    Case(cp5, cp4, [Face, Face], "intersecting", None),
+    Case(cp5, cp4, [Face, Face, Edge, Edge], "intersecting", None, True),
+    Case(
+        cp1,
+        Compound(children=cp1.get_type(Vertex)),
+        [Vertex, Vertex],
+        "mixed child type",
+        None,
+    ),
+    Case(
+        cp4,
+        Compound(children=cp3.get_type(Face)),
+        [Face, Face],
+        "mixed child type",
+        None,
+    ),
     Case(cp2, [cp3, cp4], [Edge, Edge], "multi to_intersect, intersecting", None),
-
-    Case(cv1, cp3, [Edge, Edge], "intersecting", "duplicate edges, BRepAlgoAPI_Common and _Section both return edge"),
+    Case(cv1, cp3, [Edge, Edge, Edge, Edge], "intersecting", None),  # xfail removed
     Case(sk1, cp3, [Face, Face], "intersecting", None),
     Case(pt1, cp3, [Face, Face], "intersecting", None),
-
+    Case(pt1, cp6, [Face, Face], "intersecting", None),
+    Case(pt1, cp6, [Face, Face, Edge, Edge], "intersecting", None, True),
 ]
 
-@pytest.mark.parametrize("obj, target, expected", make_params(shape_compound_matrix))
-def test_shape_compound(obj, target, expected):
-    run_test(obj, target, expected)
+@pytest.mark.parametrize("obj, target, expected, include_touched", make_params(shape_compound_matrix))
+def test_shape_compound(obj, target, expected, include_touched):
+    run_test(obj, target, expected, include_touched)
 
 # FreeCAD issue example
 c1 = CenterArc((0, 0), 10, 0, 360).edge()
@@ -437,9 +540,9 @@ freecad_matrix = [
     Case(c2, vert, [Vertex], "circle, vert, intersect", None),
 ]
 
-@pytest.mark.parametrize("obj, target, expected", make_params(freecad_matrix))
-def test_freecad(obj, target, expected):
-    run_test(obj, target, expected)
+@pytest.mark.parametrize("obj, target, expected, include_touched", make_params(freecad_matrix))
+def test_freecad(obj, target, expected, include_touched):
+    run_test(obj, target, expected, include_touched)
 
 
 # Issue tests
@@ -460,9 +563,9 @@ issues_matrix = [
     Case(e1, f1, [Edge], "issue #697", None),
 ]
 
-@pytest.mark.parametrize("obj, target, expected", make_params(issues_matrix))
-def test_issues(obj, target, expected):
-    run_test(obj, target, expected)
+@pytest.mark.parametrize("obj, target, expected, include_touched", make_params(issues_matrix))
+def test_issues(obj, target, expected, include_touched):
+    run_test(obj, target, expected, include_touched)
 
 
 # Exceptions
@@ -489,4 +592,384 @@ def make_exception_params(matrix):
 @pytest.mark.parametrize("obj, target, expected", make_exception_params(exception_matrix))
 def test_exceptions(obj, target, expected):
     with pytest.raises(Exception):
-        obj.intersect(target)
+        obj.intersect(target)
+
+
+# Direct touch() method tests
+class TestTouchMethod:
+    """Tests for direct touch() method calls to cover specific code paths."""
+
+    def test_solid_vertex_touch_on_face(self):
+        """Solid.touch(Vertex) where vertex is on a face of the solid."""
+        solid = Box(2, 2, 2)  # Box from -1 to 1 in all axes
+        # Vertex on the top face (z=1)
+        vertex = Vertex(0, 0, 1)
+        result = solid.touch(vertex)
+        assert len(result) == 1
+        assert isinstance(result[0], Vertex)
+
+    def test_solid_vertex_touch_on_edge(self):
+        """Solid.touch(Vertex) where vertex is on an edge of the solid."""
+        solid = Box(2, 2, 2)
+        # Vertex on an edge (corner of top face)
+        vertex = Vertex(1, 0, 1)
+        result = solid.touch(vertex)
+        assert len(result) == 1
+        assert isinstance(result[0], Vertex)
+
+    def test_solid_vertex_touch_on_corner(self):
+        """Solid.touch(Vertex) where vertex is on a corner of the solid."""
+        solid = Box(2, 2, 2)
+        # Vertex on a corner
+        vertex = Vertex(1, 1, 1)
+        result = solid.touch(vertex)
+        assert len(result) == 1
+        assert isinstance(result[0], Vertex)
+
+    def test_solid_vertex_touch_not_touching(self):
+        """Solid.touch(Vertex) where vertex is not on the solid boundary."""
+        solid = Box(2, 2, 2)
+        vertex = Vertex(5, 5, 5)  # Far away
+        result = solid.touch(vertex)
+        assert len(result) == 0
+
+    def test_solid_vertex_touch_inside(self):
+        """Solid.touch(Vertex) where vertex is inside the solid (not touch)."""
+        solid = Box(2, 2, 2)
+        vertex = Vertex(0, 0, 0)  # Center of box
+        result = solid.touch(vertex)
+        # Inside is not a touch - touch is boundary contact only
+        assert len(result) == 0
+
+    def test_shell_tangent_touch(self):
+        """Shell.touch(Face) for tangent contact (sphere touching plane)."""
+        # Create a hemisphere shell
+        sphere = Sphere(1).faces()[0]
+        shell = Shell([sphere])
+
+        # Create a plane tangent to the sphere at bottom (z=-1)
+        tangent_face = Face(Plane.XY.offset(-1))
+
+        result = shell.touch(tangent_face)
+        # Should find tangent vertex contact at (0, 0, -1)
+        assert len(result) >= 1
+        # Result should be vertex (tangent point)
+        vertices = [r for r in result if isinstance(r, Vertex)]
+        assert len(vertices) >= 1
+
+    def test_solid_solid_touch_faces_equal(self):
+        """Solid.touch(Solid) exercises faces_equal for duplicate face detection."""
+        b1 = Box(1, 1, 1, align=Align.MIN)
+        b2 = (
+            Box(2, 2, 0.5, align=Align.MIN)
+            - Box(1, 1.2, 1, align=Align.MIN)
+            + Pos(1, 0, 0) * Box(1, 1, 1, align=Align.MIN)
+            + Box(1, 2, 0.5, align=Align.MIN)
+        )
+        result = b1.touch(b2)
+        # Should find face contact
+        assert len(result) >= 1
+        faces = [r for r in result if isinstance(r, Face)]
+        assert len(faces) >= 1
+
+
+# ShapeList.expand() tests
+class TestShapeListExpand:
+    """Tests for ShapeList.expand() method."""
+
+    def test_expand_with_vector(self):
+        """ShapeList containing Vector objects."""
+        from build123d import Vector, ShapeList
+
+        v1 = Vector(1, 2, 3)
+        v2 = Vector(4, 5, 6)
+        shapes = ShapeList([v1, v2])
+        expanded = shapes.expand()
+        assert len(expanded) == 2
+        assert v1 in expanded
+        assert v2 in expanded
+
+    def test_expand_nested_compound(self):
+        """ShapeList with nested compounds."""
+        # Create inner compound
+        inner = Compound([Box(1, 1, 1), Pos(3, 0, 0) * Box(1, 1, 1)])
+        # Create outer compound containing inner compound
+        outer = Compound([inner, Pos(0, 3, 0) * Box(1, 1, 1)])
+
+        shapes = ShapeList([outer])
+        expanded = shapes.expand()
+
+        # Should have 3 solids after expanding nested compounds
+        solids = [s for s in expanded if isinstance(s, Solid)]
+        assert len(solids) == 3
+
+    def test_expand_shell_to_faces(self):
+        """ShapeList with Shell expands to faces."""
+        shells = Box(1, 1, 1).shells()  # Get shell from solid
+        if shells:
+            shell = shells[0]
+            shapes = ShapeList([shell])
+            expanded = shapes.expand()
+            faces = [s for s in expanded if isinstance(s, Face)]
+            assert len(faces) == 6  # Box has 6 faces
+
+    def test_expand_wire_to_edges(self):
+        """ShapeList with Wire expands to edges."""
+        wire = Rectangle(2, 2).wire()
+        shapes = ShapeList([wire])
+        expanded = shapes.expand()
+        edges = [s for s in expanded if isinstance(s, Edge)]
+        assert len(edges) == 4  # Rectangle has 4 edges
+
+    def test_expand_mixed(self):
+        """ShapeList with mixed types."""
+        from build123d import Vector
+
+        v = Vector(1, 2, 3)
+        wire = Rectangle(2, 2).wire()
+        solid = Box(1, 1, 1)
+        compound = Compound([Pos(5, 0, 0) * Box(1, 1, 1)])
+
+        shapes = ShapeList([v, wire, solid, compound])
+        expanded = shapes.expand()
+
+        # Vector stays as vector
+        assert v in expanded
+        # Wire expands to 4 edges
+        edges = [s for s in expanded if isinstance(s, Edge)]
+        assert len(edges) == 4
+        # Solid stays as solid
+        solids = [s for s in expanded if isinstance(s, Solid)]
+        assert len(solids) == 2  # Original + from compound
+
+
+class TestShellTangentTouchCoverage:
+    """Tests for Shell tangent touch to cover two_d.py lines 467-491.
+
+    These tests specifically target the Shell-specific code paths in Face.touch()
+    where we need to find which face in a Shell contains the contact point.
+    """
+
+    def test_shell_self_tangent_touch_multiple_faces(self):
+        """Shell.touch(Face) where Shell has multiple faces.
+
+        Finding face containing contact point in self Shell.
+        """
+        # Create a shell with multiple faces (half-sphere has curved + flat faces)
+        half_sphere = Sphere(1) & Pos(0, 0, 0.5) * Box(3, 3, 2)
+        shell = Shell(half_sphere.faces())
+
+        # Create a plane tangent to the curved part at x=1
+        tangent_face = Pos(1, 0, 0) * (Rot(0, 90, 0) * Rectangle(2, 2).face())
+
+        result = shell.touch(tangent_face)
+        # Should find tangent vertex at (1, 0, 0)
+        assert len(result) >= 1
+        vertices = [r for r in result if isinstance(r, Vertex)]
+        assert len(vertices) >= 1
+
+    def test_face_shell_other_tangent_touch_multiple_faces(self):
+        """Face.touch(Shell) where Shell (other) has multiple faces.
+
+        Finding face containing contact point in other Shell.
+        """
+        # Create a face
+        face = Pos(1, 0, 0) * (Rot(0, 90, 0) * Rectangle(2, 2).face())
+
+        # Create a shell with multiple faces (half-sphere)
+        half_sphere = Sphere(1) & Pos(0, 0, 0.5) * Box(3, 3, 2)
+        shell = Shell(half_sphere.faces())
+
+        result = face.touch(shell)
+        # Should find tangent vertex at (1, 0, 0)
+        assert len(result) >= 1
+        vertices = [r for r in result if isinstance(r, Vertex)]
+        assert len(vertices) >= 1
+
+    def test_shell_shell_tangent_touch_multiple_faces(self):
+        """Shell.touch(Shell) where both Shells have multiple faces.
+        """
+        # Create two half-spheres touching at their curved surfaces
+        half_sphere1 = Sphere(1) & Pos(0, 0, 0.5) * Box(3, 3, 2)
+        shell1 = Shell(half_sphere1.faces())
+
+        half_sphere2 = Pos(2, 0, 0) * (Sphere(1) & Pos(0, 0, 0.5) * Box(3, 3, 2))
+        shell2 = Shell(half_sphere2.faces())
+
+        result = shell1.touch(shell2)
+        # Should find tangent vertex at (1, 0, 0)
+        assert len(result) >= 1
+        vertices = [r for r in result if isinstance(r, Vertex)]
+        assert len(vertices) >= 1
+
+    def test_interior_tangent_contact_shell_face(self):
+        """Shell.touch(Face) with interior tangent contact (not on any edges).
+
+        Full interior tangent detection code path including Shell face
+        lookup and normal direction validation.
+
+        Contact point must be:
+        - NOT on any edge of the shell (self)
+        - NOT on any edge of the face (other)
+        """
+        import math
+
+        # Create a sphere shell
+        sphere = Sphere(2)
+        shell = Shell(sphere.faces())
+
+        # Contact at (1, 1, sqrt(2)) - away from the y=0 seam plane of the sphere
+        # This point is in the interior of the spherical surface
+        x, y, z = 1.0, 1.0, math.sqrt(2)
+
+        # Normal direction at this point on the sphere
+        normal = Vector(x, y, z).normalized()
+
+        # Create a small face tangent to sphere at this point
+        # The face must be small enough that its edges don't reach the contact point
+        tangent_plane = Plane(origin=(x, y, z), z_dir=(normal.X, normal.Y, normal.Z))
+        small_face = tangent_plane * Rectangle(0.1, 0.1).face()
+
+        result = shell.touch(small_face)
+        # Should find interior tangent vertex near (1, 1, sqrt(2))
+        assert len(result) >= 1
+        vertices = [r for r in result if isinstance(r, Vertex)]
+        assert len(vertices) >= 1
+
+    def test_interior_tangent_contact_face_shell(self):
+        """Face.touch(Shell) with interior tangent contact.
+
+        Same as above but with arguments swapped to test the 'other is Shell' path.
+        """
+        import math
+
+        # Create a sphere shell
+        sphere = Sphere(2)
+        shell = Shell(sphere.faces())
+
+        # Contact at (1, 1, sqrt(2))
+        x, y, z = 1.0, 1.0, math.sqrt(2)
+        normal = Vector(x, y, z).normalized()
+
+        # Create a small face tangent to sphere
+        tangent_plane = Plane(origin=(x, y, z), z_dir=(normal.X, normal.Y, normal.Z))
+        small_face = tangent_plane * Rectangle(0.1, 0.1).face()
+
+        # Call face.touch(shell) - 'other' is the Shell
+        result = small_face.touch(shell)
+        assert len(result) >= 1
+        vertices = [r for r in result if isinstance(r, Vertex)]
+        assert len(vertices) >= 1
+
+
+class TestSolidEdgeTangentTouch:
+    """Tests for Solid.touch(Edge) tangent cases to cover three_d.py lines 891-906.
+
+    These tests cover the BRepExtrema tangent detection for edges tangent to
+    solid surfaces (not penetrating).
+    """
+
+    def test_edge_tangent_to_cylinder(self):
+        """Edge tangent to cylinder surface returns touch vertex.
+
+        Tangent contact detection via BRepExtrema.
+        """
+        # Create a cylinder along Z axis
+        cylinder = Cylinder(1, 2)
+
+        # Create an edge that is tangent to the cylinder at x=1
+        # Edge runs along Y at x=1, z=1 (tangent to cylinder surface)
+        tangent_edge = Edge.make_line((1, -2, 1), (1, 2, 1))
+
+        result = cylinder.touch(tangent_edge)
+        # Should find tangent vertices where edge touches cylinder
+        # The edge at x=1 is tangent to the cylinder at radius=1
+        vertices = [r for r in result if isinstance(r, Vertex)]
+        # Should have at least one tangent contact point
+        assert len(vertices) >= 1
+
+    def test_edge_tangent_to_sphere(self):
+        """Edge tangent to sphere surface returns touch vertex.
+
+        Another test with spherical geometry.
+        """
+        # Create a sphere centered at origin
+        sphere = Sphere(1)
+
+        # Create an edge that is tangent to the sphere at x=1
+        # Edge runs along Z at x=1, y=0
+        tangent_edge = Edge.make_line((1, 0, -2), (1, 0, 2))
+
+        result = sphere.touch(tangent_edge)
+        # Should find tangent vertex at (1, 0, 0)
+        vertices = [r for r in result if isinstance(r, Vertex)]
+        assert len(vertices) >= 1
+
+
+class TestConvertToShapes:
+    """Tests for helpers.convert_to_shapes() to cover helpers.py."""
+
+    def test_vector_intersection(self):
+        """Shape.intersect(Vector) converts Vector to Vertex."""
+        box = Box(2, 2, 2)
+        # Vector inside the box
+        result = box.intersect(Vector(0, 0, 0))
+        assert result is not None
+        assert len(result) == 1
+        assert isinstance(result[0], Vertex)
+
+    def test_location_intersection(self):
+        """Shape.intersect(Location) converts Location to Vertex."""
+        box = Box(2, 2, 2)
+        # Location inside the box
+        result = box.intersect(Location((0, 0, 0)))
+        assert result is not None
+        assert len(result) == 1
+        assert isinstance(result[0], Vertex)
+
+    def test_location_intersection_with_rotation(self):
+        """Shape.intersect(Location with rotation) still uses position only."""
+        box = Box(2, 2, 2)
+        # Location with rotation - position is still at origin
+        loc = Location((0, 0, 0), (45, 45, 45))
+        result = box.intersect(loc)
+        assert result is not None
+        assert len(result) == 1
+        assert isinstance(result[0], Vertex)
+
+
+class TestEmptyCompoundIntersect:
+    """Tests for Compound._intersect() edge cases to cover composite.py line 741."""
+
+    def test_empty_compound_intersect(self):
+        """Empty Compound.intersect() returns None.
+
+        Early return when compound has no elements.
+        """
+        from OCP.TopoDS import TopoDS_Compound
+        from OCP.BRep import BRep_Builder
+
+        # Create an actual empty OCCT compound (has wrapped but no children)
+        builder = BRep_Builder()
+        empty_occt = TopoDS_Compound()
+        builder.MakeCompound(empty_occt)
+        empty = Compound(empty_occt)
+
+        box = Box(2, 2, 2)
+        result = empty.intersect(box)
+        assert result is None
+
+    def test_empty_compound_intersect_with_face(self):
+        """Empty Compound.intersect(Face) returns None."""
+        from OCP.TopoDS import TopoDS_Compound
+        from OCP.BRep import BRep_Builder
+
+        # Create an actual empty OCCT compound
+        builder = BRep_Builder()
+        empty_occt = TopoDS_Compound()
+        builder.MakeCompound(empty_occt)
+        empty = Compound(empty_occt)
+
+        face = Rectangle(2, 2).face()
+        result = empty.intersect(face)
+        assert result is None

tests/test_direct_api/test_location.py

@@ -405,14 +405,14 @@ class TestLocation(unittest.TestCase):
 
         self.assertIsNone(b.intersect(b.moved(Pos(X=10))))
 
-        # Look for common vertices
+        # Look for common vertices (endpoint-endpoint contacts are "touch", not "intersect")
         e1 = Edge.make_line((0, 0), (1, 0))
         e2 = Edge.make_line((1, 0), (1, 1))
         e3 = Edge.make_line((1, 0), (2, 0))
-        i = e1.intersect(e2)
+        i = e1.intersect(e2, include_touched=True)
         self.assertEqual(len(i.vertices()), 1)
         self.assertEqual(tuple(i.vertex()), (1, 0, 0))
-        i = e1.intersect(e3)
+        i = e1.intersect(e3, include_touched=True)
         self.assertEqual(len(i.vertices()), 1)
         self.assertEqual(tuple(i.vertex()), (1, 0, 0))