CSG: Rounded Cube Drilled Along Three Axes

A classic constructive solid geometry workflow expressed as a single boolean_chain call:

1. Rounded cube = cube ∩ sphere — the sphere clips the cube's corners.
2. Drilled rounded cube = result − cyl_x − cyl_y − cyl_z, where each cylinder is centered at the origin and oriented along one world axis.

All five input meshes are sent to C++ in a single boolean_chain call. Intermediate Surface_meshes never leave C++; only the final (V, F) is returned to Python.

The cylinders meet at the origin — a degenerate configuration that would trip CGAL's "Non-handled triple intersection" precondition during a naive sequential subtraction with the inexact-constructions kernel. Between corefinement steps boolean_chain applies CGAL 6.1's autorefine_triangle_soup with apply_iterative_snap_rounding(true), introduced specifically for this problem, which makes the chain robust at radii up to ~0.7 in this geometry.

from compas.geometry import Box
from compas.geometry import Cylinder
from compas.geometry import Frame
from compas.geometry import Polyhedron
from compas.geometry import Sphere
from compas_viewer import Viewer

from compas_cgal.booleans import boolean_chain


# =============================================================================
# Classic CSG: rounded cube drilled along three axes.
#   step 1: rounded cube  =  cube  ∩  sphere
#   step 2: result        =  rounded cube  −  cyl_x  −  cyl_y  −  cyl_z
#
# All five meshes are sent to C++ in a single boolean_chain() call. Each
# step's intermediate Surface_mesh stays in C++ — there is no Python ↔ C++
# round-tripping between operations. The chain runs in CGAL's exact-
# constructions kernel (EPECK), which handles the degenerate triple
# intersection at the origin where the three axis-aligned cylinders meet —
# no shifts needed. Pass `hybrid=True` to switch to the hybrid kernel
# scheme (EPICK mesh + EPECK vertex_point_map) from CGAL's "consecutive
# boolean operations with exact point maps" example.
# =============================================================================

cube = Box(2).to_vertices_and_faces(triangulated=True)
sphere = Sphere(1.3, point=[0, 0, 0]).to_vertices_and_faces(u=64, v=64, triangulated=True)


def cylinder_along(axis, radius=0.8):
    """A cylinder of length 4 oriented along a world axis.

    `Cylinder` defines its long axis as the frame's z-axis, so we build a
    frame whose z-axis matches the chosen world axis.
    """
    if axis == "x":
        frame = Frame([0, 0, 0], [0, 1, 0], [0, 0, 1])
    elif axis == "y":
        frame = Frame([0, 0, 0], [0, 0, 1], [1, 0, 0])
    else:  # "z"
        frame = Frame([0, 0, 0], [1, 0, 0], [0, 1, 0])
    return Cylinder(radius, 4.0, frame).to_vertices_and_faces(u=48, triangulated=True)


V, F = boolean_chain(
    [cube, sphere, cylinder_along("x"), cylinder_along("y"), cylinder_along("z")],
    ["intersection", "difference", "difference", "difference"],
)
shape = Polyhedron(V.tolist(), F.tolist()).to_mesh()

# =============================================================================
# Visualize
# =============================================================================

viewer = Viewer()
viewer.scene.add(shape, lineswidth=1, show_points=False)
viewer.show()