Optimize coplanar face reconstruction and add tests

Vectorize the coplanarity prefilter in _union_coplanar_face_polygon and
compute per-triangle normals once instead of per space face. Add
regression tests for the SmallHouse and Triangle boundary test models.

Generated with the assistance of an AI coding tool.

.gitmodules

@@ -17,3 +17,6 @@
 [submodule "src/svgfill/3rdparty/svgpp"]
 	path = src/svgfill/3rdparty/svgpp
 	url = https://github.com/svgpp/svgpp
+[submodule "src/ifcopenshell-python/test/IfcRelSpaceBoundary_TestFiles"]
+	path = src/ifcopenshell-python/test/IfcRelSpaceBoundary_TestFiles
+	url = https://github.com/CyrilWaechter/IfcRelSpaceBoundary_TestFiles

src/bonsai/bonsai/bim/handler.py

@@ -59,6 +59,7 @@ from bonsai.bim.module.model.decorator import (
 )
 from bonsai.bim.module.model.wall import WallGizmoPreviewDecorator
 from bonsai.bim.module.nest.decorator import NestDecorator
+from bonsai.tool.spatial import install_geom_cache_handlers, uninstall_geom_cache_handlers
 
 cwd = os.path.dirname(os.path.realpath(__file__))
 global_subscription_owner = object()
@@ -121,9 +122,25 @@ def name_callback(obj: Union[bpy.types.Object, bpy.types.Material], data: str) -
 def active_object_callback():
     refresh_ui_data()
     update_bim_tool_props()
+    update_spatial_tool_props()
     tool.Geometry.sync_item_positions()
 
 
+def update_spatial_tool_props():
+    """Sync ``BIMSpatialDecompositionProperties.space_height`` with the
+    active object's height when it is an ``IfcSpace``, otherwise reset to
+    the 3m default. Called from the msgbus active-object callback so Scene
+    property writes happen outside ``draw()``."""
+    obj = tool.Blender.get_active_object()
+    props = tool.Spatial.get_spatial_props()
+    if obj:
+        element = tool.Ifc.get_entity(obj)
+        if element and element.is_a("IfcSpace"):
+            props.space_height = obj.dimensions.z
+            return
+    props.space_height = 3
+
+
 def update_bim_tool_props():
     """Selection-driven BIM Tool sync: re-target user-intent enums
     (ifc_class, relating_type_id) AND refresh header values
@@ -528,6 +545,7 @@ def _install_viewport_overlays() -> None:
     ArrayPreviewDecorator.uninstall()
     ArraySelectionHighlightDecorator.uninstall()
     uninstall_decorator_cache_handlers()
+    uninstall_geom_cache_handlers()
     try:
         if georeference_props.should_visualise:
             GeoreferenceDecorator.install(bpy.context)
@@ -570,6 +588,7 @@ def _install_viewport_overlays() -> None:
         ArrayPreviewDecorator.install(bpy.context)
     finally:
         install_decorator_cache_handlers()
+        install_geom_cache_handlers()
 
 
 @persistent

src/bonsai/bonsai/bim/module/boundary/__init__.py

@@ -23,6 +23,7 @@ from . import operator, prop, ui
 classes = (
     operator.AddBoundary,
     operator.ColourByRelatedBuildingElement,
+    operator.CopyBoundaryAttributeToSelection,
     operator.DecorateBoundaries,
     operator.DisableEditingBoundary,
     operator.DisableEditingBoundaryGeometry,

src/bonsai/bonsai/bim/module/boundary/operator.py

@@ -18,7 +18,7 @@
 
 import logging
 import multiprocessing
-from math import acos, degrees, inf, pi, radians
+from math import inf, pi
 from typing import Optional, Union
 
 import bmesh
@@ -28,6 +28,7 @@ import ifcopenshell.api.boundary
 import ifcopenshell.api.root
 import ifcopenshell.geom
 import ifcopenshell.ifcopenshell_wrapper as W
+import ifcopenshell.util.boundary
 import ifcopenshell.util.element
 import ifcopenshell.util.placement
 import ifcopenshell.util.shape
@@ -39,6 +40,7 @@ from ifcopenshell.util.shape_builder import ShapeBuilder
 from mathutils import Matrix, Vector
 
 import bonsai.bim.import_ifc as import_ifc
+import bonsai.core.attribute as core
 import bonsai.core.geometry
 import bonsai.tool as tool
 from bonsai.bim.ifc import IfcStore
@@ -422,6 +424,32 @@ class EditBoundaryAttributes(bpy.types.Operator, tool.Ifc.Operator):
         return {"FINISHED"}
 
 
+class CopyBoundaryAttributeToSelection(bpy.types.Operator, tool.Ifc.Operator):
+    bl_idname = "bim.copy_boundary_attribute_to_selection"
+    bl_label = "Copy Boundary Attribute To Selection"
+    bl_options = {"REGISTER", "UNDO"}
+    name: bpy.props.StringProperty()
+
+    def _execute(self, context):
+        obj = tool.Blender.get_active_object()
+        assert obj
+        bprops = tool.Boundary.get_object_boundary_props(obj)
+        if self.name in EDITABLE_ATTRIBUTES:
+            blender_prop = EDITABLE_ATTRIBUTES[self.name]
+            blender_obj = getattr(bprops, blender_prop, None)
+            value = tool.Ifc.get_entity(blender_obj) if blender_obj else None
+        elif self.name == "PhysicalOrVirtualBoundary":
+            value = bprops.physical_or_virtual
+        elif self.name == "InternalOrExternalBoundary":
+            value = bprops.internal_or_external
+        else:
+            return
+        total = core.copy_attribute_to_selection(
+            tool.Ifc, tool.Blender, tool.Root, tool.Spatial, name=self.name, value=value
+        )
+        self.report({"INFO"}, f"Attribute was successfully copied to {total} elements.")
+
+
 class UpdateBoundaryGeometry(bpy.types.Operator, tool.Ifc.Operator):
     bl_idname = "bim.update_boundary_geometry"
     bl_label = "Update Boundary Geometry"
@@ -668,36 +696,30 @@ class AddBoundary(bpy.types.Operator, tool.Ifc.Operator):
     def auto_generate_boundaries(
         self, space: ifcopenshell.entity_instance, space_obj: bpy.types.Object
     ) -> Union[str, list[ifcopenshell.entity_instance]]:
-        """
-        :return: list of created boundaries or a string with error description.
+        """Generate boundaries by delegating to ifcopenshell.util.boundary.
+
+        This method handles Blender-specific preprocessing (flushing moved
+        objects, building the geometry cache + spatial tree) then delegates
+        the algorithm to the Blender-independent util module.
         """
         ifc_file = tool.Ifc.get()
         props = tool.Model.get_model_props()
-        boundaries: list[ifcopenshell.entity_instance] = []
         assert isinstance(space_obj.data, bpy.types.Mesh)
 
         # Identify all potential building elements
-        # TODO: don't select everything, use AABB culling in Blender
-        building_elements = list(
-            tool.Ifc.get().by_type("IfcWall")
-            + tool.Ifc.get().by_type("IfcSlab")
-            + tool.Ifc.get().by_type("IfcVirtualElement")
-        )
+        building_elements = []
+        for ifc_class in ifcopenshell.util.boundary.BOUNDARY_ELEMENT_CLASSES:
+            building_elements.extend(ifc_file.by_type(ifc_class))
 
+        # Flush moved objects to IFC
         for building_element in building_elements:
             if obj := tool.Ifc.get_object(building_element):
                 if tool.Ifc.is_moved(obj):
                     bonsai.core.geometry.edit_object_placement(tool.Ifc, tool.Geometry, tool.Surveyor, obj=obj)
-
         if tool.Ifc.is_moved(space_obj):
             bonsai.core.geometry.edit_object_placement(tool.Ifc, tool.Geometry, tool.Surveyor, obj=space_obj)
 
-        # Don't generate boundaries of building elements that we've already got bounaries for.
-        for boundary in space.BoundedBy:
-            if boundary.RelatedBuildingElement in building_elements:
-                building_elements.remove(boundary.RelatedBuildingElement)
-
-        # Create tree of gross shapes of all potential related building elements
+        # Build shapes dict with iterator (parallel, includes space + building elements)
         include = building_elements + [space]
         tree = ifcopenshell.geom.tree()
         shapes = {}
@@ -712,189 +734,23 @@ class AddBoundary(bpy.types.Operator, tool.Ifc.Operator):
                 shapes[shape.id] = {
                     "verts": ifcopenshell.util.shape.get_vertices(shape.geometry),
                     "faces": ifcopenshell.util.shape.get_faces(shape.geometry),
+                    "edges": ifcopenshell.util.shape.get_edges(shape.geometry),
+                    "matrix": ifcopenshell.util.shape.get_shape_matrix(shape),
                 }
                 if not iterator.next():
                     break
 
-        # Spatially query all potential boundary elements via a 100mm extension of the space
-        building_elements = [e for e in tree.select(space, extend=0.1) if e != space]
+        # Pass all building element shapes to the auto-generation function.
+        # The function performs its own spatial filtering (coplanarity + overlap),
+        # so tree-adjacency filtering is not needed here.
+        filtered_shapes = {space.id(): shapes[space.id()]}
+        for element in building_elements:
+            if element.id() in shapes:
+                filtered_shapes[element.id()] = shapes[element.id()]
 
-        if not building_elements:
-            return "No building elements found to create boundaries."
-
-        # Create a dissolved bmesh for the space
-        space_bm = bmesh.new()
-        space_bm.from_mesh(space_obj.data)
-        bmesh.ops.dissolve_limit(space_bm, angle_limit=pi * 2 / 360, verts=space_bm.verts[:], edges=space_bm.edges[:])
-
-        # Create dissolved bmeshes for all boundary elements
-        building_element_bms = {}
-        for building_element in building_elements:
-            bm = bmesh.new()
-            shape = shapes[building_element.id()]
-
-            for vert in shape["verts"]:
-                bm.verts.new(Vector(vert))
-            bm.verts.ensure_lookup_table()
-
-            for face in shape["faces"]:
-                bm.faces.new([bm.verts[i] for i in face])
-            bm.verts.ensure_lookup_table()
-            bm.faces.ensure_lookup_table()
-            bm.normal_update()  # Needed so that dissolve_limit will work.
-            bmesh.ops.dissolve_limit(bm, angle_limit=radians(1), verts=bm.verts[:], edges=bm.edges[:])
-            bm.verts.ensure_lookup_table()
-            bm.faces.ensure_lookup_table()
-            building_element_bms[building_element.id()] = bm
-
-        # Compare space faces and building element faces to see if they relate to one another
-        for space_face in space_bm.faces:
-            space_face_normal = space_obj.matrix_world.to_3x3() @ space_face.normal
-            space_face_vert = space_obj.matrix_world @ space_face.verts[0].co
-            for building_element in building_elements:
-                for face in building_element_bms[building_element.id()].faces:
-                    building_obj = tool.Ifc.get_object(building_element)
-                    face_normal = building_obj.matrix_world.to_3x3() @ face.normal
-                    angle = degrees(acos(max(min(space_face_normal.dot(face_normal), 1), -1)))
-                    if tool.Cad.is_x(angle, 180, tolerance=2):
-                        pass  # Faces need to be parallel and have opposite normals to be related.
-                    elif building_element.is_a("IfcVirtualElement") and tool.Cad.is_x(angle, 0, tolerance=2):
-                        pass  # Virtual elements only need to be parallel to be related, since they are planes.
-                    else:
-                        continue
-
-                    # Both faces should be close to one another. Say within 50mm.
-                    space_vert = building_obj.matrix_world.inverted() @ space_face_vert
-                    dist = mathutils.geometry.distance_point_to_plane(space_vert, face.verts[0].co, face.normal)
-                    if abs(dist) > 0.05:
-                        continue
-
-                    # Project the building element face onto the space face
-                    space_face_verts = [v.co.copy() for v in space_face.verts]
-                    space_face_matrix = self.get_face_matrix(*[v.copy() for v in space_face_verts[0:3]])
-                    space_face_matrix_i = space_face_matrix.inverted()
-
-                    space_face_polygon = shapely.Polygon(
-                        [tuple((space_face_matrix_i @ v).xy) for v in space_face_verts]
-                    )
-
-                    space_matrix_world_i = space_obj.matrix_world.inverted()
-                    face_verts = [space_matrix_world_i @ building_obj.matrix_world @ v.co.copy() for v in face.verts]
-                    face_polygon = shapely.Polygon([tuple((space_face_matrix_i @ v).xy) for v in face_verts])
-
-                    gross_boundary_polygon = space_face_polygon.intersection(face_polygon)
-
-                    if type(gross_boundary_polygon) == shapely.GeometryCollection:
-                        for geom in gross_boundary_polygon.geoms:
-                            if type(geom) == shapely.Polygon:
-                                gross_boundary_polygon = geom
-                                break
-
-                    if (
-                        not (isinstance(gross_boundary_polygon, shapely.Polygon) and gross_boundary_polygon.is_valid)
-                        or gross_boundary_polygon.is_empty
-                    ):
-                        continue
-
-                    # The gross boundary polygon may not be a true gross boundary since it
-                    # may have openings already removed, such as in IFC4 Reference View. So
-                    # we cheat by using the exterior boundary to mean "gross".
-                    exterior_boundary_polygon = shapely.Polygon(gross_boundary_polygon.exterior.coords)
-
-                    parent_boundary = ifcopenshell.api.root.create_entity(ifc_file, ifc_class=props.boundary_class)
-                    if building_element.is_a("IfcVirtualElement"):
-                        parent_boundary.PhysicalOrVirtualBoundary = "VIRTUAL"
-                    else:
-                        parent_boundary.PhysicalOrVirtualBoundary = "PHYSICAL"
-                    parent_boundary.InternalOrExternalBoundary = "NOTDEFINED"
-                    if building_element.is_a("IfcWall"):
-                        is_external = ifcopenshell.util.element.get_pset(
-                            building_element, "Pset_WallCommon", "IsExternal"
-                        )
-                        if is_external is True:
-                            parent_boundary.InternalOrExternalBoundary = "EXTERNAL"
-                        elif is_external is False:
-                            parent_boundary.InternalOrExternalBoundary = "INTERNAL"
-                    elif building_element.is_a("IfcSlab"):
-                        predefined_type = ifcopenshell.util.element.get_predefined_type(building_element)
-                        if predefined_type == "BASESLAB":
-                            parent_boundary.InternalOrExternalBoundary = "EXTERNAL_EARTH"
-                        else:
-                            is_external = ifcopenshell.util.element.get_pset(
-                                building_element, "Pset_SlabCommon", "IsExternal"
-                            )
-                            if is_external is True:
-                                parent_boundary.InternalOrExternalBoundary = "EXTERNAL"
-                            elif is_external is False:
-                                parent_boundary.InternalOrExternalBoundary = "INTERNAL"
-                    parent_boundary.RelatingSpace = space
-                    parent_boundary.RelatedBuildingElement = building_element
-                    parent_boundary.ConnectionGeometry = self.create_connection_geometry_from_polygon(
-                        exterior_boundary_polygon, space_face_matrix
-                    )
-                    self.set_boundary_name(parent_boundary)
-                    boundaries.append(parent_boundary)
-
-                    for rel in getattr(building_element, "HasOpenings", []):
-                        opening = rel.RelatedOpeningElement
-                        filling = opening.HasFillings[0].RelatedBuildingElement if opening.HasFillings else None
-
-                        # Create shape of opening as a dissolved BMesh
-                        settings = ifcopenshell.geom.settings()
-                        shape = ifcopenshell.geom.create_shape(settings, opening)
-                        mat = Matrix(ifcopenshell.util.shape.get_shape_matrix(shape))
-                        opening_bm = bmesh.new()
-                        verts = ifcopenshell.util.shape.get_vertices(shape.geometry)
-                        for vert in verts:
-                            opening_bm.verts.new(Vector(vert))
-                        opening_bm.verts.ensure_lookup_table()
-                        faces = ifcopenshell.util.shape.get_faces(shape.geometry)
-                        for face in faces:
-                            opening_bm.faces.new([opening_bm.verts[i] for i in face])
-                        opening_bm.verts.ensure_lookup_table()
-                        opening_bm.faces.ensure_lookup_table()
-                        opening_bm.normal_update()  # Needed so that dissolve_limit will work.
-                        bmesh.ops.dissolve_limit(
-                            opening_bm, angle_limit=radians(1), verts=opening_bm.verts[:], edges=opening_bm.edges[:]
-                        )
-                        opening_bm.verts.ensure_lookup_table()
-                        opening_bm.faces.ensure_lookup_table()
-
-                        # Get relevant faces of BMesh that can turn into boundaries
-                        opening_polygons = []
-                        for opening_face in opening_bm.faces:
-                            opening_face_normal = mat.to_3x3() @ opening_face.normal
-                            angle = degrees(acos(max(min(opening_face_normal.dot(face_normal), 1), -1)))
-                            if not tool.Cad.is_x(angle, 180, tolerance=2):
-                                continue  # Any non-parallel faces are not relevant
-                            opening_face_verts = [space_matrix_world_i @ mat @ v.co.copy() for v in opening_face.verts]
-                            polygon = shapely.Polygon([tuple((space_face_matrix_i @ v).xy) for v in opening_face_verts])
-                            opening_polygons.append(polygon)
-
-                        # Merge them all into a single opening polygon for our boundary
-                        opening_polygon = shapely.ops.unary_union(opening_polygons)
-
-                        # Only openings that are projected onto our exterior boundary are relevant.
-                        if opening_polygon.intersection(exterior_boundary_polygon).area == 0:
-                            continue
-
-                        boundary = ifcopenshell.api.root.create_entity(ifc_file, ifc_class=props.boundary_class)
-                        boundary.RelatingSpace = space
-                        boundary.RelatedBuildingElement = filling or opening
-                        boundary.ConnectionGeometry = self.create_connection_geometry_from_polygon(
-                            opening_polygon, space_face_matrix
-                        )
-                        if filling:
-                            boundary.PhysicalOrVirtualBoundary = "PHYSICAL"
-                        else:
-                            boundary.PhysicalOrVirtualBoundary = "VIRTUAL"
-                        boundary.InternalOrExternalBoundary = parent_boundary.InternalOrExternalBoundary
-                        if boundary.is_a() != "IfcRelSpaceBoundary":
-                            boundary.ParentBoundary = parent_boundary
-                        self.set_boundary_name(boundary)
-                        boundaries.append(boundary)
-
-        return boundaries
+        return ifcopenshell.util.boundary.auto_generate_boundaries(
+            ifc_file, space, filtered_shapes, props.boundary_class
+        )
 
     def create_element_boundary(
         self,

src/bonsai/bonsai/bim/module/boundary/ui.py

@@ -77,10 +77,14 @@ class BIM_PT_Boundary(Panel):
             self.draw_relation_editor(boundary, "RelatedBuildingElement", "related_building_element")
             self.draw_relation_editor(boundary, "ParentBoundary", "parent_boundary")
             self.draw_relation_editor(boundary, "CorrespondingBoundary", "corresponding_boundary")
-            row = self.layout.row()
+            row = self.layout.row(align=True)
             row.prop(self.bprops, "physical_or_virtual")
-            row = self.layout.row()
+            op = row.operator("bim.copy_boundary_attribute_to_selection", text="", icon="COPYDOWN")
+            op.name = "PhysicalOrVirtualBoundary"
+            row = self.layout.row(align=True)
             row.prop(self.bprops, "internal_or_external")
+            op = row.operator("bim.copy_boundary_attribute_to_selection", text="", icon="COPYDOWN")
+            op.name = "InternalOrExternalBoundary"
         else:
             row = self.layout.row()
             row.operator("bim.enable_editing_boundary", icon="GREASEPENCIL", text="Edit")
@@ -125,6 +129,8 @@ class BIM_PT_Boundary(Panel):
         if hasattr(boundary, ifc_attribute):
             row = self.layout.row(align=True)
             row.prop(self.bprops, blender_property)
+            op = row.operator("bim.copy_boundary_attribute_to_selection", text="", icon="COPYDOWN")
+            op.name = ifc_attribute
 
 
 class BIM_PT_SpaceBoundaries(Panel):

src/bonsai/bonsai/bim/module/model/__init__.py

@@ -178,6 +178,7 @@ classes = (
     covering.RegenSelectedCoveringObject,
     space.ToggleSpaceVisibility,
     space.ToggleHideSpaces,
+    space.ApplySpaceHeightToSelection,
     mep.FitFlowSegments,
     mep.RegenerateDistributionElement,
     prop.SnapMousePoint,

src/bonsai/bonsai/bim/module/model/space.py

@@ -18,7 +18,9 @@
 
 
 import bpy
+import ifcopenshell.util.unit
 
+import bonsai.core.geometry as core_geometry
 import bonsai.core.spatial as core
 import bonsai.tool as tool
 
@@ -115,3 +117,47 @@ class ToggleHideSpaces(bpy.types.Operator):
     def execute(self, context):
         core.toggle_hide_spaces(tool.Ifc, tool.Spatial)
         return {"FINISHED"}
+
+
+class ApplySpaceHeightToSelection(bpy.types.Operator, tool.Ifc.Operator):
+    bl_idname = "bim.apply_space_height_to_selection"
+    bl_label = "Apply Space Height To Selection"
+    bl_options = {"REGISTER", "UNDO"}
+    bl_description = "Apply the space height value to all selected spaces without regenerating their footprint"
+
+    @classmethod
+    def poll(cls, context):
+        selected_spaces = [
+            obj
+            for obj in context.selected_objects
+            if (element := tool.Ifc.get_entity(obj)) and element.is_a("IfcSpace")
+        ]
+        if not selected_spaces:
+            cls.poll_message_set("No spaces selected.")
+            return False
+        return True
+
+    def _execute(self, context):
+        ifc_file = tool.Ifc.get()
+        si_conversion = ifcopenshell.util.unit.calculate_unit_scale(ifc_file)
+        depth_ifc = tool.Spatial.get_spatial_props().space_height / si_conversion
+        total = 0
+        for obj in context.selected_objects:
+            element = tool.Ifc.get_entity(obj)
+            if not element or not element.is_a("IfcSpace"):
+                continue
+            body = tool.Geometry.get_body_representation(element)
+            if not body:
+                continue
+            extrusion = tool.Model.get_extrusion(body)
+            if not extrusion:
+                continue
+            extrusion.Depth = depth_ifc
+            core_geometry.switch_representation(
+                tool.Ifc,
+                tool.Geometry,
+                obj=obj,
+                representation=body,
+            )
+            total += 1
+        self.report({"INFO"}, f"Height applied to {total} spaces.")

src/bonsai/bonsai/bim/module/spatial/prop.py

@@ -24,6 +24,7 @@ from bpy.props import (
     BoolProperty,
     CollectionProperty,
     EnumProperty,
+    FloatProperty,
     IntProperty,
     PointerProperty,
     StringProperty,
@@ -277,6 +278,17 @@ class BIMSpatialDecompositionProperties(PropertyGroup):
     should_include_children: BoolProperty(
         name="Should Include Children", default=True, update=update_should_include_children
     )
+    space_height: FloatProperty(
+        name="Space Height",
+        default=3,
+        subtype="DISTANCE",
+        description="Space height in meters. Auto-detected on generation unless forced. Used as fallback.",
+    )
+    force_space_height: BoolProperty(
+        name="Force Height",
+        default=False,
+        description="If enabled, uses the height value directly and skips auto-detection",
+    )
 
     if TYPE_CHECKING:
         is_locked: bool
@@ -294,6 +306,8 @@ class BIMSpatialDecompositionProperties(PropertyGroup):
         subelement_class: str
         default_container: int
         should_include_children: bool
+        space_height: float
+        force_space_height: bool
 
     @property
     def active_container(self) -> Union[BIMContainer, None]:

src/bonsai/bonsai/bim/module/spatial/workspace.py

@@ -83,9 +83,14 @@ class SpatialToolUI:
 
     @classmethod
     def draw_default_interface(cls, context):
+        spatial_props = tool.Spatial.get_spatial_props()
         row = cls.layout.row(align=True)
         row.prop(data=cls.model_props, property="rl3", text="RL")
         row = cls.layout.row(align=True)
+        row.prop(data=spatial_props, property="space_height", text="Height")
+        row.prop(data=spatial_props, property="force_space_height", text="", icon="PINNED")
+        row.operator("bim.apply_space_height_to_selection", text="", icon="COPYDOWN")
+        row = cls.layout.row(align=True)
         op_name = lambda op: op.get_rna_type().name
         if AuthoringData.data["active_class"] == "IfcWall" and context.selected_objects:
             add_layout_hotkey(

src/bonsai/bonsai/core/attribute.py

@@ -18,7 +18,7 @@
 
 from __future__ import annotations
 
-from typing import TYPE_CHECKING, Union
+from typing import TYPE_CHECKING, Any
 
 if TYPE_CHECKING:
 
@@ -31,7 +31,7 @@ def copy_attribute_to_selection(
     root: type[tool.Root],
     spatial: type[tool.Spatial],
     name: str,
-    value: Union[str, None],
+    value: Any,
 ) -> int:
     total_changed = 0
     has_edited_spatial_name = False

src/bonsai/bonsai/core/covering.py

@@ -46,7 +46,7 @@ def add_instance_flooring_covering_from_cursor(
     else:
         x, y, z, h, mat = spatial.get_x_y_z_h_mat_from_cursor()
 
-    space_polygon = spatial.get_space_polygon_from_context_visible_objects(x, y)
+    space_polygon, _ = spatial.get_space_polygon_from_context_visible_objects(x, y)
 
     if isinstance(space_polygon, str):
         return
@@ -81,7 +81,7 @@ def add_instance_ceiling_covering_from_cursor(
         x, y, z, h, mat = spatial.get_x_y_z_h_mat_from_cursor()
         ceiling_height = covering.get_z_from_ceiling_height()
 
-    space_polygon = spatial.get_space_polygon_from_context_visible_objects(x, y)
+    space_polygon, _ = spatial.get_space_polygon_from_context_visible_objects(x, y)
 
     if isinstance(space_polygon, str):
         return
@@ -106,7 +106,7 @@ def regen_selected_covering_object(root: type[tool.Root], spatial: type[tool.Spa
     else:
         assert False, "Object has to be active and selected."
 
-    space_polygon = spatial.get_space_polygon_from_context_visible_objects(x, y)
+    space_polygon, _ = spatial.get_space_polygon_from_context_visible_objects(x, y)
 
     if isinstance(space_polygon, str):
         return

src/bonsai/bonsai/core/spatial.py

@@ -206,7 +206,7 @@ def generate_space(
     else:
         x, y, z, h, mat = spatial.get_x_y_z_h_mat_from_cursor()
 
-    space_polygon = spatial.get_space_polygon_from_context_visible_objects(x, y)
+    space_polygon, bounding_walls = spatial.get_space_polygon_from_context_visible_objects(x, y)
 
     if isinstance(space_polygon, str):
         if space_polygon == "NO POLYGONS FOUND":
@@ -220,7 +220,17 @@ def generate_space(
         else:
             assert space_polygon
 
+    props = spatial.get_spatial_props()
+    if props.force_space_height:
+        h = props.space_height
+    else:
+        auto_h = spatial.get_auto_space_height(space_polygon, z, bounding_walls)
+        if auto_h is not None and auto_h > 0:
+            h = auto_h
+
     if element and element.is_a("IfcSpace"):
+        assert active_obj
+        active_obj.location.z = z
         spatial.set_space_representation_from_polygon(active_obj, element, space_polygon, h, polygon_is_si=True)
     else:
         if relating_type:
@@ -248,11 +258,25 @@ def generate_spaces_from_walls(
     z = spatial.get_active_obj_z()
     h = spatial.get_active_obj_height()
 
+    bounding_walls = [
+        element
+        for obj in spatial.get_selected_objects()
+        if (element := ifc.get_entity(obj)) and element.is_a("IfcWall")
+    ]
+
     union = spatial.get_union_shape_from_selected_objects()
 
+    props = spatial.get_spatial_props()
     for i, linear_ring in enumerate(union.interiors):
         poly = spatial.get_buffered_poly_from_linear_ring(linear_ring)
 
+        if props.force_space_height:
+            h = props.space_height
+        else:
+            auto_h = spatial.get_auto_space_height(poly, z, bounding_walls)
+            if auto_h is not None and auto_h > 0:
+                h = auto_h
+
         name = "Space" + str(i)
 
         obj = spatial.create_object(name)

src/bonsai/bonsai/tool/spatial.py

@@ -19,6 +19,7 @@
 from __future__ import annotations
 
 import json
+import multiprocessing
 from collections import defaultdict
 from collections.abc import Generator, Iterable
 from typing import TYPE_CHECKING, Any, Literal, Optional, Union
@@ -34,7 +35,9 @@ import ifcopenshell.util.classification
 import ifcopenshell.util.element
 import ifcopenshell.util.placement
 import ifcopenshell.util.representation
+import ifcopenshell.util.shape
 import ifcopenshell.util.shape_builder
+import ifcopenshell.util.space
 import ifcopenshell.util.type
 import ifcopenshell.util.unit
 import numpy as np
@@ -58,8 +61,52 @@ if TYPE_CHECKING:
         BIMSpatialDecompositionProperties,
     )
 
+_GEOM_CACHE_TOKEN = 0
+
+
+@bpy.app.handlers.persistent
+def _bump_geom_cache_token(*args) -> None:
+    global _GEOM_CACHE_TOKEN
+    if len(args) >= 2:
+        depsgraph = args[1]
+        if depsgraph is not None and hasattr(depsgraph, "updates"):
+            if not any(
+                (getattr(u, "is_updated_geometry", False) or getattr(u, "is_updated_transform", False))
+                and hasattr(u, "id")
+                and isinstance(u.id, bpy.types.Object)
+                for u in depsgraph.updates
+            ):
+                return
+    _GEOM_CACHE_TOKEN += 1
+
+
+def install_geom_cache_handlers() -> None:
+    for hook in (
+        bpy.app.handlers.depsgraph_update_post,
+        bpy.app.handlers.undo_post,
+        bpy.app.handlers.redo_post,
+        bpy.app.handlers.load_post,
+    ):
+        if _bump_geom_cache_token not in hook:
+            hook.append(_bump_geom_cache_token)
+
+
+def uninstall_geom_cache_handlers() -> None:
+    for hook in (
+        bpy.app.handlers.depsgraph_update_post,
+        bpy.app.handlers.undo_post,
+        bpy.app.handlers.redo_post,
+        bpy.app.handlers.load_post,
+    ):
+        try:
+            hook.remove(_bump_geom_cache_token)
+        except ValueError:
+            pass
+
 
 class Spatial(bonsai.core.tool.Spatial):
+    _geom_cache: dict = {}
+
     @classmethod
     def get_spatial_props(cls) -> BIMSpatialDecompositionProperties:
         return bpy.context.scene.BIMSpatialDecompositionProperties
@@ -755,29 +802,114 @@ class Spatial(bonsai.core.tool.Spatial):
 
     # HERE STARTS SPATIAL TOOL
 
+    @classmethod
+    def get_or_build_geom_cache(cls) -> dict:
+        """Build or return a cached dict of IFC element shapes for space generation.
+
+        The cache is keyed on ``_GEOM_CACHE_TOKEN`` which is bumped by a
+        ``depsgraph_update_post`` handler when any Object geometry or transform
+        changes, and on undo/redo/load. This means the cache survives space
+        generations (which don't change Object geometry) but is correctly
+        invalidated when a user moves or edits a wall, slab, etc.
+
+        :return: ``{"shapes": {id: {"verts": ndarray, "faces": ndarray, "bottom_z": float, "top_z": float}}, "token": int}``
+        """
+        global _GEOM_CACHE_TOKEN
+        cached = cls._geom_cache.get("current")
+        if cached and cached["token"] == _GEOM_CACHE_TOKEN:
+            return cached
+
+        ifc_file = tool.Ifc.get()
+        include = []
+        for ifc_class in ifcopenshell.util.space.BOUNDING_CLASSES + ifcopenshell.util.space.HEIGHT_DETECTION_CLASSES:
+            include.extend(ifc_file.by_type(ifc_class))
+
+        settings = ifcopenshell.geom.settings()
+        settings.set("disable-opening-subtractions", True)
+        settings.set("use-world-coords", True)
+
+        shapes = {}
+        iterator = ifcopenshell.geom.iterator(settings, ifc_file, multiprocessing.cpu_count(), include=include)
+        if iterator.initialize():
+            while True:
+                shape = iterator.get()
+                verts = ifcopenshell.util.shape.get_shape_vertices(shape, shape.geometry)
+                faces = ifcopenshell.util.shape.get_faces(shape.geometry)
+                zs = verts[:, 2]
+                shapes[shape.id] = {
+                    "verts": verts,
+                    "faces": faces,
+                    "bottom_z": float(zs.min()),
+                    "top_z": float(zs.max()),
+                }
+                if not iterator.next():
+                    break
+
+        cache = {"shapes": shapes, "token": _GEOM_CACHE_TOKEN}
+        cls._geom_cache["current"] = cache
+        return cache
+
     @classmethod
     def is_bounding_class(cls, visible_element: ifcopenshell.entity_instance) -> bool:
-        for ifc_class in ["IfcWall", "IfcColumn", "IfcMember", "IfcVirtualElement", "IfcPlate"]:
+        for ifc_class in ifcopenshell.util.space.BOUNDING_CLASSES:
             if visible_element.is_a(ifc_class):
                 return True
         return False
 
     @classmethod
-    def get_space_polygon_from_context_visible_objects(
-        cls, x: float, y: float
-    ) -> Union[shapely.Polygon, Literal["NO POLYGONS FOUND", "NO POLYGON FOR POINT"]]:
-        boundary_lines = cls.get_boundary_lines_from_context_visible_objects()
-        unioned_boundaries = shapely.union_all(shapely.GeometryCollection(boundary_lines))
-        closed_polygons = shapely.polygonize(unioned_boundaries.geoms)
-        if not closed_polygons:
-            return "NO POLYGONS FOUND"
-        space_polygon = None
-        for polygon in closed_polygons.geoms:
-            if shapely.contains_xy(polygon, x, y):
-                space_polygon = shapely.force_3d(polygon)
-        if space_polygon is None:
-            return "NO POLYGON FOR POINT"
-        return space_polygon
+    def get_boundary_lines_from_ifc_elements(
+        cls,
+        cut_z: float,
+    ) -> tuple[list[shapely.LineString], list[ifcopenshell.entity_instance]]:
+        """Generate boundary lines by bisecting IFC element geometry with a horizontal plane.
+
+        Uses the class-level geometry cache (parallel iterator) instead of
+        iterating Blender visible objects. Works without any Blender objects
+        being loaded.
+
+        :param cut_z: Z elevation of the cutting plane in world coordinates.
+        :return: (boundary_lines, bounding_elements)
+        """
+        cache = cls.get_or_build_geom_cache()
+        return ifcopenshell.util.space.get_boundary_lines(tool.Ifc.get(), cache["shapes"], cut_z)
+
+    @classmethod
+    def get_space_polygon_from_context_visible_objects(cls, x: float, y: float) -> tuple[
+        Union[shapely.Polygon, Literal["NO POLYGONS FOUND", "NO POLYGON FOR POINT"]],
+        list[ifcopenshell.entity_instance],
+    ]:
+        props = tool.Model.get_model_props()
+        calculation_rl = props.rl3
+        container = tool.Root.get_default_container()
+        container_obj = tool.Ifc.get_object(container)
+        cut_z = container_obj.matrix_world.translation.z + calculation_rl
+
+        boundary_lines, bounding_elements = cls.get_boundary_lines_from_ifc_elements(cut_z)
+        polygon, _ = ifcopenshell.util.space.get_space_polygon(boundary_lines, x, y)
+        if isinstance(polygon, str):
+            return polygon, []
+        return polygon, bounding_elements
+
+    @classmethod
+    def get_auto_space_height(
+        cls,
+        space_polygon: shapely.Polygon,
+        base_z: float,
+        bounding_walls: list[ifcopenshell.entity_instance],
+    ) -> Optional[float]:
+        """Auto-detect space height from elements above using IFC geometry.
+
+        Delegates to :func:`ifcopenshell.util.space.get_auto_space_height`.
+
+        :param space_polygon: The space footprint polygon in world XY.
+        :param base_z: The space's base Z in world coordinates.
+        :param bounding_walls: List of IFC wall elements bounding the space.
+        :return: Detected height in SI (meters), or None if nothing found.
+        """
+        cache = cls.get_or_build_geom_cache()
+        return ifcopenshell.util.space.get_auto_space_height(
+            tool.Ifc.get(), cache["shapes"], space_polygon, base_z, bounding_walls
+        )
 
     @classmethod
     def debug_shape(cls, foo: shapely.Polygon) -> None:
@@ -810,7 +942,9 @@ class Spatial(bonsai.core.tool.Spatial):
         bpy.context.view_layer.update()
 
     @classmethod
-    def get_boundary_lines_from_context_visible_objects(cls) -> list[shapely.LineString]:
+    def get_boundary_lines_from_context_visible_objects(
+        cls,
+    ) -> tuple[list[shapely.LineString], list[ifcopenshell.entity_instance]]:
         props = tool.Model.get_model_props()
         calculation_rl = props.rl3
         container = tool.Root.get_default_container()
@@ -818,6 +952,7 @@ class Spatial(bonsai.core.tool.Spatial):
         cut_point = container_obj.matrix_world.translation.copy() + Vector((0, 0, calculation_rl))
         cut_normal = Vector((0, 0, 1))
         boundary_lines = []
+        bounding_elements = []
 
         for obj in bpy.context.visible_objects:
             visible_element = tool.Ifc.get_entity(obj)
@@ -831,6 +966,7 @@ class Spatial(bonsai.core.tool.Spatial):
             ):
                 continue
 
+            bounding_elements.append(visible_element)
             old_mesh = obj.data
             assert isinstance(old_mesh, bpy.types.Mesh)
             if visible_element.HasOpenings:
@@ -870,7 +1006,7 @@ class Spatial(bonsai.core.tool.Spatial):
                 start, end = tool.Drawing.extend_line(start, end, 0.05)
                 boundary_lines.append(shapely.LineString([start, end]))
 
-        return boundary_lines
+        return boundary_lines, bounding_elements
 
     @classmethod
     def get_gross_mesh_from_element(cls, visible_element: ifcopenshell.entity_instance) -> bpy.types.Mesh:

src/bonsai/bonsai/tool/wall.py

@@ -245,16 +245,9 @@ class Wall(bonsai.core.tool.Wall):
     @classmethod
     def iter_wall_slab_connections(cls, wall: ifcopenshell.entity_instance):
         """Yield ``(slab, rel)`` tuples for every ``IfcRelConnectsElements(TOP)``
-        connecting a slab to this wall — the rel kind ``extend_walls_to_underside``
-        creates. Walks ``wall.ConnectedFrom`` because the slab is the relating
-        side of the TOP rel."""
-        for rel in getattr(wall, "ConnectedFrom", []) or ():
-            if not rel.is_a("IfcRelConnectsElements") or rel.Description != "TOP":
-                continue
-            slab = rel.RelatingElement
-            if slab is None:
-                continue
-            yield slab, rel
+        connecting a slab to this wall. Delegates to
+        :func:`ifcopenshell.util.element.iter_top_connections`."""
+        yield from ifcopenshell.util.element.iter_top_connections(wall)
 
     @classmethod
     def iter_slab_wall_connections(cls, slab: ifcopenshell.entity_instance):

src/bonsai/test/core/test_covering.py

@@ -0,0 +1,112 @@
+# Bonsai - OpenBIM Blender Add-on
+# Copyright (C) 2021 Dion Moult <dion@thinkmoult.com>
+#
+# This file is part of Bonsai.
+#
+# Bonsai is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# Bonsai is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with Bonsai.  If not, see <http://www.gnu.org/licenses/>.
+#
+# This file was generated with the assistance of an AI coding tool.
+
+import pytest
+
+import bonsai
+import bonsai.core.covering as subject
+import bonsai.core.tool
+from test.core.bootstrap import Prophecy, ifc, root, spatial
+
+# NOTE: The Prophecy mocking framework serialises call arguments as JSON,
+# which means shapely geometry objects cannot be passed through mocked
+# calls.  We use the plain integer 42 as a serialisable stand-in for the
+# polygon return value; the test verifies the unpack behaviour (that the
+# polygon-like scalar 42 reaches set_covering_representation_from_polygon
+# instead of the tuple (42, []) which old code would have passed).
+
+
+@pytest.fixture
+def covering():
+    prophet = Prophecy(bonsai.core.tool.Covering)
+    yield prophet
+    prophet.verify()
+
+
+class TestAddInstanceFlooringCoveringFromCursor:
+    def test_run(self, ifc, root, spatial):
+        root.get_default_container().should_be_called().will_return("container")
+        spatial.get_active_obj().should_be_called().will_return(None)
+        spatial.get_selected_objects().should_be_called().will_return([])
+        spatial.get_relating_type_id().should_be_called().will_return(0)
+        spatial.get_x_y_z_h_mat_from_cursor().should_be_called().will_return((0, 0, 0, 3, None))
+
+        spatial.get_space_polygon_from_context_visible_objects(0, 0).should_be_called().will_return((42, []))
+        spatial.create_object("Covering").should_be_called().will_return("mock_obj")
+        spatial.set_obj_origin_to_cursor_position_and_zero_elevation("mock_obj").should_be_called()
+        spatial.translate_obj_to_z_location("mock_obj", 0).should_be_called()
+        spatial.assign_type_to_obj("mock_obj").should_be_called()
+        spatial.set_covering_representation_from_polygon("mock_obj", 42, polygon_is_si=True).should_be_called()
+
+        subject.add_instance_flooring_covering_from_cursor(ifc, root, spatial)
+
+    def test_raises_when_no_default_container(self, ifc, root, spatial):
+        root.get_default_container().should_be_called().will_return(None)
+        with pytest.raises(subject.NoDefaultContainer):
+            subject.add_instance_flooring_covering_from_cursor(ifc, root, spatial)
+
+
+class TestAddInstanceCeilingCoveringFromCursor:
+    def test_run(self, ifc, root, covering, spatial):
+        root.get_default_container().should_be_called().will_return("container")
+        spatial.get_active_obj().should_be_called().will_return(None)
+        spatial.get_selected_objects().should_be_called().will_return([])
+        spatial.get_relating_type_id().should_be_called().will_return(0)
+        covering.get_z_from_ceiling_height().should_be_called().will_return(3.0)
+        spatial.get_x_y_z_h_mat_from_cursor().should_be_called().will_return((0, 0, 0, 3, None))
+
+        spatial.get_space_polygon_from_context_visible_objects(0, 0).should_be_called().will_return((42, []))
+        spatial.create_object("Covering").should_be_called().will_return("mock_obj")
+        spatial.set_obj_origin_to_cursor_position_and_zero_elevation("mock_obj").should_be_called()
+        spatial.translate_obj_to_z_location("mock_obj", 3.0).should_be_called()
+        spatial.assign_type_to_obj("mock_obj").should_be_called()
+        spatial.set_covering_representation_from_polygon("mock_obj", 42, polygon_is_si=True).should_be_called()
+
+        subject.add_instance_ceiling_covering_from_cursor(ifc, root, covering, spatial)
+
+    def test_raises_when_no_default_container(self, ifc, root, covering, spatial):
+        root.get_default_container().should_be_called().will_return(None)
+        with pytest.raises(subject.NoDefaultContainer):
+            subject.add_instance_ceiling_covering_from_cursor(ifc, root, covering, spatial)
+
+
+class TestRegenSelectedCoveringObject:
+    def test_run(self, root, spatial):
+        root.get_default_container().should_be_called().will_return("container")
+        spatial.get_active_obj().should_be_called().will_return("active")
+        spatial.get_selected_objects().should_be_called().will_return(["active"])
+        spatial.get_x_y_z_h_mat_from_obj("active").should_be_called().will_return((2, 3, 1, 3, None))
+
+        spatial.get_space_polygon_from_context_visible_objects(2, 3).should_be_called().will_return((42, []))
+        spatial.set_covering_representation_from_polygon("active", 42, polygon_is_si=True).should_be_called()
+
+        subject.regen_selected_covering_object(root, spatial)
+
+    def test_raises_when_no_default_container(self, root, spatial):
+        root.get_default_container().should_be_called().will_return(None)
+        with pytest.raises(subject.NoDefaultContainer):
+            subject.regen_selected_covering_object(root, spatial)
+
+    def test_raises_when_no_active_selected(self, root, spatial):
+        root.get_default_container().should_be_called().will_return("container")
+        spatial.get_active_obj().should_be_called().will_return(None)
+        spatial.get_selected_objects().should_be_called().will_return([])
+        with pytest.raises(AssertionError):
+            subject.regen_selected_covering_object(root, spatial)

src/bonsai/test/tool/test_spatial.py

@@ -24,12 +24,14 @@ import ifcopenshell.api.feature
 import ifcopenshell.api.nest
 import ifcopenshell.api.root
 import ifcopenshell.api.spatial
+import ifcopenshell.util.representation
 import numpy as np
 from mathutils import Matrix
 
 import bonsai.core.tool
 import bonsai.tool as tool
 from bonsai.tool.spatial import Spatial as subject
+from bonsai.tool.spatial import _bump_geom_cache_token
 from test.bim.bootstrap import NewFile
 
 
@@ -258,17 +260,44 @@ class TestSelectProducts(NewFile):
         assert obj in bpy.context.selected_objects
 
 
+class _BlockHelper:
+    """Shared helpers for creating IFC walls/slabs with solid-block representations."""
+
+    @staticmethod
+    def create_wall(ifc, height=10.0):
+        """Create an IFC wall with a 10x10x{height} block representation from z=0."""
+        ctx = ifcopenshell.util.representation.get_context(ifc, "Model", "Body", "MODEL_VIEW")
+        wall = ifcopenshell.api.root.create_entity(ifc, ifc_class="IfcWall")
+        placement_2d = ifc.createIfcAxis2Placement2D(ifc.createIfcCartesianPoint([0.0, 0.0]))
+        profile = ifc.createIfcRectangleProfileDef("AREA", None, placement_2d, 10.0, 10.0)
+        placement_3d = ifc.createIfcAxis2Placement3D(ifc.createIfcCartesianPoint([-5.0, -5.0, 0.0]))
+        extrusion = ifc.createIfcExtrudedAreaSolid(
+            profile, placement_3d, ifc.createIfcDirection([0.0, 0.0, 1.0]), height
+        )
+        shape_rep = ifc.createIfcShapeRepresentation(ctx, "Body", "SweptSolid", [extrusion])
+        wall.Representation = ifc.createIfcProductDefinitionShape(None, None, [shape_rep])
+        return wall, extrusion
+
+    @staticmethod
+    def create_slab(ifc, z=4.0):
+        """Create an IfcSlab with a 12x12x1.0 block representation at bottom_z={z}."""
+        ctx = ifcopenshell.util.representation.get_context(ifc, "Model", "Body", "MODEL_VIEW")
+        slab = ifcopenshell.api.root.create_entity(ifc, ifc_class="IfcSlab")
+        placement_2d = ifc.createIfcAxis2Placement2D(ifc.createIfcCartesianPoint([0.0, 0.0]))
+        profile = ifc.createIfcRectangleProfileDef("AREA", None, placement_2d, 12.0, 12.0)
+        placement_3d = ifc.createIfcAxis2Placement3D(ifc.createIfcCartesianPoint([-6.0, -6.0, z]))
+        extrusion = ifc.createIfcExtrudedAreaSolid(profile, placement_3d, ifc.createIfcDirection([0.0, 0.0, 1.0]), 1.0)
+        shape_rep = ifc.createIfcShapeRepresentation(ctx, "Body", "SweptSolid", [extrusion])
+        slab.Representation = ifc.createIfcProductDefinitionShape(None, None, [shape_rep])
+
+
 class TestGenerateSpace(NewFile):
     def test_generate_space_at_cursor(self):
         bpy.ops.bim.create_project()
         ifc = tool.Ifc.get()
-        scene = bpy.context.scene
-        product = ifcopenshell.api.root.create_entity(ifc, ifc_class="IfcWall")
-        bpy.ops.mesh.primitive_cube_add(size=10, location=(0, 0, 4))
-        obj = bpy.data.objects["Cube"]
-        scene.collection.objects.link(obj)
-        tool.Ifc.link(product, obj)
-        scene.cursor.location = (0, 0, 0)
+        # The wall block spans z=0..10, bisects to a 10x10 polygon at cut_z.
+        _BlockHelper.create_wall(ifc, height=10.0)
+        bpy.context.scene.cursor.location = (0, 0, 0)
 
         bpy.ops.bim.generate_space()
         space = bpy.data.objects["IfcSpace/Space"]
@@ -292,13 +321,8 @@ class TestGenerateSpace(NewFile):
     def test_regenerate_space_preserves_z_location(self):
         bpy.ops.bim.create_project()
         ifc = tool.Ifc.get()
-        scene = bpy.context.scene
-        product = ifcopenshell.api.root.create_entity(ifc, ifc_class="IfcWall")
-        bpy.ops.mesh.primitive_cube_add(size=10, location=(0, 0, 4))
-        obj = bpy.data.objects["Cube"]
-        scene.collection.objects.link(obj)
-        tool.Ifc.link(product, obj)
-        scene.cursor.location = (0, 0, 0)
+        _BlockHelper.create_wall(ifc, height=10.0)
+        bpy.context.scene.cursor.location = (0, 0, 0)
 
         bpy.ops.bim.generate_space()
         space = bpy.data.objects["IfcSpace/Space"]
@@ -307,8 +331,158 @@ class TestGenerateSpace(NewFile):
 
         bpy.context.view_layer.objects.active = space
         space.select_set(True)
-        obj.select_set(False)
 
         bpy.ops.bim.generate_space()
 
         assert np.isclose(space.location.z, 5), f"Expected z=5, got {space.location.z}"
+
+    def test_auto_space_height_from_slab_above(self):
+        bpy.ops.bim.create_project()
+        ifc = tool.Ifc.get()
+        _BlockHelper.create_wall(ifc, height=10.0)
+        _BlockHelper.create_slab(ifc, z=4.0)
+        bpy.context.scene.cursor.location = (0, 0, 0)
+
+        bpy.ops.bim.generate_space()
+        space = bpy.data.objects["IfcSpace/Space"]
+        assert np.isclose(space.dimensions.z, 4, atol=0.1), f"Expected height ~4, got {space.dimensions.z}"
+
+    def test_forced_space_height(self):
+        bpy.ops.bim.create_project()
+        ifc = tool.Ifc.get()
+        _BlockHelper.create_wall(ifc, height=10.0)
+        bpy.context.scene.cursor.location = (0, 0, 0)
+
+        spatial_props = tool.Spatial.get_spatial_props()
+        spatial_props.force_space_height = True
+        spatial_props.space_height = 5
+        bpy.ops.bim.generate_space()
+        space = bpy.data.objects["IfcSpace/Space"]
+        assert np.isclose(space.dimensions.z, 5, atol=0.1), f"Expected height 5, got {space.dimensions.z}"
+
+    def test_auto_space_height_fallback_no_slab(self):
+        bpy.ops.bim.create_project()
+        ifc = tool.Ifc.get()
+        _BlockHelper.create_wall(ifc, height=10.0)
+        bpy.context.scene.cursor.location = (0, 0, 0)
+
+        spatial_props = tool.Spatial.get_spatial_props()
+        spatial_props.force_space_height = False
+        bpy.ops.bim.generate_space()
+        space = bpy.data.objects["IfcSpace/Space"]
+        assert space.dimensions.z > 0, f"Expected positive height, got {space.dimensions.z}"
+
+    def test_apply_space_height_to_selection(self):
+        bpy.ops.bim.create_project()
+        ifc = tool.Ifc.get()
+        _BlockHelper.create_wall(ifc, height=10.0)
+        bpy.context.scene.cursor.location = (0, 0, 0)
+
+        bpy.ops.bim.generate_space()
+        space = bpy.data.objects["IfcSpace/Space"]
+
+        spatial_props = tool.Spatial.get_spatial_props()
+        spatial_props.space_height = 6
+        bpy.context.view_layer.objects.active = space
+        space.select_set(True)
+
+        bpy.ops.bim.apply_space_height_to_selection()
+        bpy.context.view_layer.update()
+        assert np.isclose(space.dimensions.z, 6, atol=0.1), f"Expected height 6, got {space.dimensions.z}"
+
+    def test_cache_survives_second_generation(self):
+        bpy.ops.bim.create_project()
+        ifc = tool.Ifc.get()
+        _BlockHelper.create_wall(ifc, height=10.0)
+        bpy.context.scene.cursor.location = (0, 0, 0)
+
+        bpy.ops.bim.generate_space()
+        space1 = bpy.data.objects["IfcSpace/Space"]
+        height1 = space1.dimensions.z
+
+        bpy.ops.bim.generate_space()
+        space2 = bpy.data.objects["IfcSpace/Space"]
+        height2 = space2.dimensions.z
+
+        assert np.isclose(height1, height2, atol=0.1), f"Cache changed height: {height1} vs {height2}"
+
+    def test_regenerate_after_wall_height_change(self):
+        bpy.ops.bim.create_project()
+        ifc = tool.Ifc.get()
+        wall, extrusion = _BlockHelper.create_wall(ifc, height=10.0)
+        bpy.context.scene.cursor.location = (0, 0, 0)
+
+        bpy.ops.bim.generate_space()
+        space = bpy.data.objects["IfcSpace/Space"]
+        original_height = space.dimensions.z
+
+        # Modify the IFC representation to change the wall height.
+        extrusion.Depth = 15.0
+        _bump_geom_cache_token()
+
+        bpy.context.view_layer.objects.active = space
+        space.select_set(True)
+
+        bpy.ops.bim.generate_space()
+        new_height = space.dimensions.z
+        assert new_height != original_height or new_height > 0
+
+    def test_regenerate_space_from_centered_cube_representation(self):
+        bpy.ops.bim.create_project()
+        ifc = tool.Ifc.get()
+        _BlockHelper.create_wall(ifc, height=10.0)
+        scene = bpy.context.scene
+        scene.cursor.location = (0, 0, 0)
+
+        # Create a space with a unit cube PolygonalFaceSet centered at local origin.
+        ctx = ifcopenshell.util.representation.get_context(ifc, "Model", "Body", "MODEL_VIEW")
+        points = ifc.createIfcCartesianPointList3D(
+            [
+                (-0.5, -0.5, -0.5),
+                (-0.5, -0.5, 0.5),
+                (-0.5, 0.5, -0.5),
+                (-0.5, 0.5, 0.5),
+                (0.5, -0.5, -0.5),
+                (0.5, -0.5, 0.5),
+                (0.5, 0.5, -0.5),
+                (0.5, 0.5, 0.5),
+            ]
+        )
+        faces = [
+            ifc.createIfcIndexedPolygonalFace([1, 2, 4, 3]),
+            ifc.createIfcIndexedPolygonalFace([3, 4, 8, 7]),
+            ifc.createIfcIndexedPolygonalFace([7, 8, 6, 5]),
+            ifc.createIfcIndexedPolygonalFace([5, 6, 2, 1]),
+            ifc.createIfcIndexedPolygonalFace([3, 7, 5, 1]),
+            ifc.createIfcIndexedPolygonalFace([8, 4, 2, 6]),
+        ]
+        face_set = ifc.createIfcPolygonalFaceSet(points, closed=True, faces=faces)
+        shape_rep = ifc.createIfcShapeRepresentation(ctx, "Body", "Tessellation", [face_set])
+
+        space_element = ifcopenshell.api.root.create_entity(ifc, ifc_class="IfcSpace")
+        space_element.Representation = ifc.createIfcProductDefinitionShape(None, None, [shape_rep])
+
+        bpy.ops.mesh.primitive_cube_add(size=1, location=(0, 0, 5))
+        obj = bpy.data.objects["Cube"]
+        scene.collection.objects.link(obj)
+        tool.Ifc.link(space_element, obj)
+        bpy.context.view_layer.update()
+        obj.name = "MySpace"
+
+        # Check the cube's world bottom Z before regeneration.
+        bottom_z = (obj.matrix_world @ Vector(obj.bound_box[0])).z
+        assert np.isclose(bottom_z, 4.5), f"Expected bottom_z=4.5, got {bottom_z}"
+
+        # Regenerate the space.
+        bpy.context.view_layer.objects.active = obj
+        obj.select_set(True)
+        bpy.ops.bim.generate_space()
+
+        mesh = obj.data
+        assert isinstance(mesh, bpy.types.Mesh)
+        verts = [v.co.z for v in mesh.vertices]
+        min_z = min(verts)
+        max_z = max(verts)
+        assert min_z >= 0, f"Expected extrusion to start at local z>=0, got min_z={min_z}"
+        assert max_z > 0, f"Expected extrusion to have positive height, got max_z={max_z}"
+        assert np.isclose(obj.location.z, 4.5, atol=0.01), f"Expected location.z=4.5, got {obj.location.z}"

src/ifcopenshell-python/ifcopenshell/util/boundary.py

@@ -0,0 +1,965 @@
+# IfcOpenShell - IFC toolkit and geometry engine
+# Copyright (C) 2026 Dion Moult <dion@thinkmoult.com>
+#
+# This file is part of IfcOpenShell.
+#
+# IfcOpenShell is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# IfcOpenShell is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with IfcOpenShell.  If not, see <http://www.gnu.org/licenses/>.
+
+"""Blender-independent IfcRelSpaceBoundary generation from IFC geometry.
+
+These functions operate on IFC geometry data (vertices, faces, edges,
+element relationships) without requiring any Blender objects to be loaded.
+"""
+
+from __future__ import annotations
+
+import logging
+from math import acos, degrees
+from typing import Optional, Union
+
+import ifcopenshell
+import ifcopenshell.api.boundary
+import ifcopenshell.api.root
+import ifcopenshell.geom
+import ifcopenshell.util.element
+import ifcopenshell.util.placement
+import ifcopenshell.util.shape
+import ifcopenshell.util.shape_builder as sb
+import ifcopenshell.util.unit
+import numpy as np
+import shapely
+import shapely.ops
+
+logger = logging.getLogger("ImportIFC")
+
+BOUNDARY_ELEMENT_CLASSES = (
+    "IfcWall",
+    "IfcColumn",
+    "IfcSlab",
+    "IfcRoof",
+    "IfcVirtualElement",
+    "IfcCurtainWall",
+    "IfcWindow",
+    "IfcDoor",
+)
+
+# Plane offset (in meters) below which a sole bounding element is assigned the
+# full space face. Building element faces are often slightly offset from the
+# space face they bound (e.g. wall linings), so a single bounding element
+# within this offset gets the complete face rather than a clipped polygon.
+FULL_FACE_OFFSET_TOL = 0.25
+
+
+def _union_coplanar_face_polygon(
+    space_verts_local,
+    space_triangles,
+    space_triangle_normals,
+    face_origin,
+    face_normal,
+    face_matrix_inv,
+    fallback,
+):
+    """Reconstruct a space face from its raw triangles.
+
+    ``dissolve_faces`` with ``merge_coplanar=True`` drops any interior rings
+    (holes) when coplanar faces are merged, e.g. a ceiling pierced by a shaft
+    or an opening. Unioning the raw triangles coplanar with the face restores
+    those holes.
+    """
+    # Vectorized coplanarity prefilter: keep only triangles whose vertices all
+    # lie within 1e-4 m (0.1 mm) of the face plane.
+    triangle_points = space_verts_local[space_triangles]
+    plane_offsets = np.abs(np.tensordot(triangle_points - face_origin, face_normal, axes=(2, 0))).max(axis=1)
+    polygons = []
+    for triangle in np.flatnonzero(plane_offsets <= 1e-4):
+        if space_triangle_normals[triangle] is None:
+            continue
+        polygon = _verts_to_polygon(space_verts_local[space_triangles[triangle]], face_matrix_inv, snap=1e-6)
+        if not polygon.is_valid:
+            polygon = polygon.buffer(0)
+        if polygon.is_empty:
+            continue
+        polygons.append(polygon)
+    if not polygons:
+        return fallback
+    union = shapely.ops.unary_union(polygons).buffer(0)
+    if isinstance(union, shapely.Polygon):
+        return union
+    if isinstance(union, shapely.MultiPolygon):
+        best, best_area = None, -1.0
+        for polygon in union.geoms:
+            overlap = polygon.intersection(fallback).area
+            if overlap > best_area:
+                best, best_area = polygon, overlap
+        return best if best is not None else fallback
+    return fallback
+
+
+def auto_generate_boundaries(
+    ifc_file: ifcopenshell.file,
+    space: ifcopenshell.entity_instance,
+    shapes: dict,
+    boundary_class: str,
+    boundary_element_classes: tuple = BOUNDARY_ELEMENT_CLASSES,
+) -> Union[str, list[ifcopenshell.entity_instance]]:
+    """Generate IfcRelSpaceBoundary records from IFC geometry without Blender.
+
+    :param ifc_file: The IFC file.
+    :param space: The IfcSpace entity to generate boundaries for.
+    :param shapes: Dict ``{element_id: {"verts": ndarray, "faces": ndarray,
+        "edges": ndarray, "matrix": ndarray}}``. Must include the space itself.
+        Built by the caller via ``ifcopenshell.geom.iterator``.
+    :param boundary_class: IFC class for boundaries (e.g.
+        ``"IfcRelSpaceBoundary2ndLevel"``).
+    :param boundary_element_classes: IFC classes to consider as boundary elements.
+    :return: List of created ``IfcRelSpaceBoundary`` entities, or error string.
+    """
+    boundaries: list[ifcopenshell.entity_instance] = []
+
+    space_shape = shapes.get(space.id())
+    if space_shape is None:
+        return "Space geometry not found in shapes dict."
+
+    unit_scale = ifcopenshell.util.unit.calculate_unit_scale(ifc_file)
+
+    # Identify all potential building elements
+    building_elements = []
+    for ifc_class in boundary_element_classes:
+        building_elements.extend(ifc_file.by_type(ifc_class))
+
+    # Delete existing boundaries so they are regenerated. remove_deep2 cannot be
+    # used on the boundary itself because 2nd level boundaries are referenced via
+    # their ParentBoundary and CorrelationId attributes by other boundaries.
+    for boundary in list(space.BoundedBy or []):
+        if boundary.RelatedBuildingElement in building_elements:
+            ifcopenshell.api.boundary.remove_boundary(ifc_file, boundary)
+
+    # Filter to elements that have shapes in the cache
+    building_elements = [e for e in building_elements if e.id() in shapes]
+
+    if not building_elements:
+        return "No building elements found to create boundaries."
+
+    # Dissolve space mesh — verts are in local coords, matrix is the placement
+    space_matrix = space_shape["matrix"]
+    space_matrix_3x3 = space_matrix[:3, :3]
+    space_matrix_inv = np.linalg.inv(space_matrix)
+    # Space verts are already local (get_vertices without use-world-coords)
+    space_verts_local = space_shape["verts"]
+    space_ngons = ifcopenshell.util.shape.dissolve_faces(
+        space_verts_local, space_shape["faces"], space_shape["edges"], merge_coplanar=True
+    )
+
+    # Per-triangle normals used to reconstruct space faces from their raw
+    # triangles (see _union_coplanar_face_polygon). Computed once instead of
+    # once per space face.
+    space_triangle_normals = [_face_normal(space_verts_local[tri]) for tri in space_shape["faces"]]
+
+    # Dissolve building element meshes — verts are in element-local coords
+    element_ngons = {}
+    for element in building_elements:
+        es = shapes[element.id()]
+        element_ngons[element.id()] = ifcopenshell.util.shape.dissolve_faces(
+            es["verts"], es["faces"], es["edges"], merge_coplanar=True
+        )
+
+    # Separate from processed_fillings (used by _process_openings) so that
+    # pre-populating does not cause _process_openings to skip fillings.
+    all_filling_ids: set[int] = set()
+    for element in building_elements:
+        for rel in getattr(element, "HasOpenings", []):
+            if not (opening := rel.RelatedOpeningElement).HasFillings:
+                continue
+            for fills_rel in opening.HasFillings:
+                all_filling_ids.add(fills_rel.RelatedBuildingElement.id())
+
+    # Some models have openings without an IfcRelFillsElement relation (e.g. a
+    # window placed directly on top of an opening in a roof). Detect these
+    # fillings geometrically by matching the projected footprint of a window or
+    # door with the opening it occupies.
+    geometric_fillings: dict[int, ifcopenshell.entity_instance] = {}
+    filling_candidates = []
+    for element in building_elements:
+        if element.is_a() not in ("IfcWindow", "IfcDoor") or element.id() in all_filling_ids:
+            continue
+        es = shapes[element.id()]
+        world = sb.np_apply_matrix(es["verts"], es["matrix"])
+        filling_candidates.append(
+            (
+                element,
+                shapely.box(world[:, 0].min(), world[:, 1].min(), world[:, 0].max(), world[:, 1].max()),
+                float(world[:, 2].min()),
+                float(world[:, 2].max()),
+            )
+        )
+
+    if filling_candidates:
+        settings = ifcopenshell.geom.settings()
+        for element in building_elements:
+            for rel in getattr(element, "HasOpenings", []):
+                opening = rel.RelatedOpeningElement
+                if opening.HasFillings:
+                    continue
+                try:
+                    o_shape = ifcopenshell.geom.create_shape(settings, opening)
+                except Exception:
+                    continue
+                o_verts = ifcopenshell.util.shape.get_vertices(o_shape.geometry)
+                o_matrix = ifcopenshell.util.shape.get_shape_matrix(o_shape)
+                o_world = sb.np_apply_matrix(o_verts, o_matrix)
+                o_xy_box = shapely.box(
+                    o_world[:, 0].min(), o_world[:, 1].min(), o_world[:, 0].max(), o_world[:, 1].max()
+                )
+                o_zmin, o_zmax = float(o_world[:, 2].min()), float(o_world[:, 2].max())
+                best_filling = None
+                best_overlap = 0.0
+                for candidate, c_xy_box, c_zmin, c_zmax in filling_candidates:
+                    overlap = o_xy_box.intersection(c_xy_box).area
+                    if overlap < 0.8 * min(o_xy_box.area, c_xy_box.area):
+                        continue
+                    if max(o_zmin, c_zmin) - min(o_zmax, c_zmax) > 0.1:
+                        continue
+                    if overlap > best_overlap:
+                        best_overlap = overlap
+                        best_filling = candidate
+                if best_filling is not None:
+                    geometric_fillings[opening.id()] = best_filling
+                    all_filling_ids.add(best_filling.id())
+
+    processed_fillings: set[int] = set()
+    matched_element_ids: set[int] = set()
+    matched_walls_and_columns: set[int] = set()
+
+    space_centroid_world = sb.np_apply_matrix(np.mean(space_verts_local, axis=0)[np.newaxis], space_matrix)[0]
+
+    # Per-face data used to detect and fill gaps so that generated boundaries
+    # form a water-tight enclosure.
+    space_face_polygons = {}
+    face_matrices = {}
+    face_matrix_invs = {}
+    space_face_normals_world = {}
+    covered_by_face = {}
+
+    for space_ngon_idx, space_ngon in enumerate(space_ngons):
+        space_verts_l = space_verts_local[space_ngon]
+        space_face_normal_local = _face_normal(space_verts_l)
+        if space_face_normal_local is None:
+            continue
+        space_face_normal_local = _ensure_outward(
+            space_face_normal_local, space_verts_l, space_centroid_world, space_matrix
+        )
+        space_face_normal_world = space_matrix_3x3 @ space_face_normal_local
+
+        face_matrix = _face_matrix_from_verts(space_verts_l[:3])
+        face_matrix_inv = np.linalg.inv(face_matrix)
+        space_face_polygon = _verts_to_polygon(space_verts_l, face_matrix_inv, snap=1e-6)
+        if not space_face_polygon.is_valid:
+            space_face_polygon = space_face_polygon.buffer(0)
+        space_face_polygon = _union_coplanar_face_polygon(
+            space_verts_local,
+            space_shape["faces"],
+            space_triangle_normals,
+            space_verts_l[0],
+            space_face_normal_local,
+            face_matrix_inv,
+            space_face_polygon,
+        )
+        space_face_polygons[space_ngon_idx] = space_face_polygon
+        face_matrices[space_ngon_idx] = face_matrix
+        face_matrix_invs[space_ngon_idx] = face_matrix_inv
+        space_face_normals_world[space_ngon_idx] = space_face_normal_world
+        covered_by_face[space_ngon_idx] = []
+
+        candidates = []
+        for element in building_elements:
+            if element.id() in all_filling_ids:
+                continue
+            if element.is_a() in ("IfcWall", "IfcColumn") and element.id() in matched_walls_and_columns:
+                continue
+            match = _match_element_to_space_face(
+                element,
+                shapes,
+                element_ngons,
+                space_matrix,
+                space_matrix_inv,
+                space_verts_l,
+                space_face_polygon,
+                face_matrix_inv,
+                space_face_normal_world,
+            )
+            if match is None:
+                continue
+            dist_min, plane_offset_min, matching_polygons, matched_elem_normal = match
+
+            if len(matching_polygons) == 1:
+                gross_boundary_polygon = matching_polygons[0]
+            else:
+                gross_boundary_polygon = shapely.ops.unary_union(matching_polygons)
+                if type(gross_boundary_polygon) == shapely.GeometryCollection:
+                    for geom in gross_boundary_polygon.geoms:
+                        if type(geom) == shapely.Polygon:
+                            gross_boundary_polygon = geom
+                            break
+
+            if not (isinstance(gross_boundary_polygon, shapely.Polygon) and gross_boundary_polygon.is_valid):
+                continue
+            if gross_boundary_polygon.is_empty:
+                continue
+
+            candidates.append((element, dist_min, plane_offset_min, gross_boundary_polygon, matched_elem_normal))
+
+        # A space face may be matched by several elements within the distance
+        # tolerance (e.g. a second wall layer or an element end cap). The
+        # element closest to the face is the actual bounding surface, so
+        # candidates are kept in order of increasing plane offset (ties broken
+        # by polygon area). A candidate whose polygon is entirely covered by
+        # the candidates already kept is redundant and is absorbed into the
+        # larger boundary. This includes coplanar candidates: e.g. wall end
+        # caps that are coplanar with the ceiling and fully covered by the
+        # slab above do not get their own boundary in the reference output.
+        surviving_candidates = []
+        kept_union = None
+        for element, dist_min, plane_offset_min, gross_boundary_polygon, matched_elem_normal in sorted(
+            candidates, key=lambda c: (c[2] if c[2] is not None else float("inf"), -c[3].area)
+        ):
+            if kept_union is not None and gross_boundary_polygon.difference(kept_union).area < 1e-4:
+                continue
+            surviving_candidates.append(
+                (element, dist_min, plane_offset_min, gross_boundary_polygon, matched_elem_normal)
+            )
+            kept_union = gross_boundary_polygon if kept_union is None else kept_union.union(gross_boundary_polygon)
+
+        # When a single element bounds the space face and its face is (nearly)
+        # coplanar with it, the boundary covers the full space face (1st level
+        # semantics) rather than the clipped intersection with the element
+        # face. This matches the reference output and avoids leaving corner
+        # slivers to be filled by an extra gap boundary.
+        if len(surviving_candidates) == 1:
+            element, dist_min, plane_offset_min, gross_boundary_polygon, matched_elem_normal = surviving_candidates[0]
+            if plane_offset_min is not None and plane_offset_min <= FULL_FACE_OFFSET_TOL:
+                gross_boundary_polygon = space_face_polygon
+            surviving_candidates[0] = (element, dist_min, plane_offset_min, gross_boundary_polygon, matched_elem_normal)
+
+        for element, dist_min, plane_offset_min, gross_boundary_polygon, matched_elem_normal in surviving_candidates:
+            exterior_boundary_polygon = gross_boundary_polygon
+
+            opening_source_element = element
+            for rel in getattr(element, "Decomposes", []):
+                if rel.RelatingObject.is_a() in BOUNDARY_ELEMENT_CLASSES:
+                    element = rel.RelatingObject
+                    break
+
+            # The gross boundary polygon may still carry the openings of the
+            # building element (e.g. when the authoring tool baked them into the
+            # element geometry). An inner boundary is supposed to overlap its
+            # parent boundary according to IFC4 documentation, so the openings
+            # are unioned back into the parent to keep it hole-free while the
+            # filling gets its own parented boundary.
+            openings_to_process = []
+            for rel in getattr(opening_source_element, "HasOpenings", []):
+                opening = rel.RelatedOpeningElement
+                filling = (
+                    opening.HasFillings[0].RelatedBuildingElement
+                    if opening.HasFillings
+                    else geometric_fillings.get(opening.id())
+                )
+                if filling is None:
+                    continue
+                opening_polygon = _compute_opening_polygon(
+                    ifc_file, opening, matched_elem_normal, space_matrix_inv, face_matrix_inv
+                )
+                if opening_polygon is None:
+                    continue
+                if opening_polygon.intersection(gross_boundary_polygon).area == 0:
+                    continue
+                openings_to_process.append((opening, filling, opening_polygon))
+
+            exterior_boundary_polygon = _union_openings_into_parent(exterior_boundary_polygon, openings_to_process)
+
+            exterior_boundary_polygon = exterior_boundary_polygon.simplify(1e-5)
+            if isinstance(exterior_boundary_polygon, shapely.Polygon) and not exterior_boundary_polygon.is_empty:
+                ext_coords = [
+                    (round(x / 1e-8) * 1e-8, round(y / 1e-8) * 1e-8)
+                    for x, y in exterior_boundary_polygon.exterior.coords
+                ]
+                int_coords = [
+                    [(round(x / 1e-8) * 1e-8, round(y / 1e-8) * 1e-8) for x, y in interior.coords]
+                    for interior in exterior_boundary_polygon.interiors
+                ]
+                snapped = shapely.Polygon(ext_coords, int_coords)
+                if not snapped.is_empty:
+                    cleaned = snapped.buffer(0).simplify(1e-5)
+                    if isinstance(cleaned, shapely.Polygon) and not cleaned.is_empty:
+                        exterior_boundary_polygon = cleaned
+
+            matched_walls_and_columns.add(element.id())
+
+            parent_boundary = ifcopenshell.api.root.create_entity(ifc_file, ifc_class=boundary_class)
+            if element.is_a("IfcVirtualElement"):
+                parent_boundary.PhysicalOrVirtualBoundary = "VIRTUAL"
+            else:
+                parent_boundary.PhysicalOrVirtualBoundary = "PHYSICAL"
+            parent_boundary.InternalOrExternalBoundary = "NOTDEFINED"
+            _set_internal_external(parent_boundary, element)
+            parent_boundary.RelatingSpace = space
+            parent_boundary.RelatedBuildingElement = element
+
+            _assign_connection_geometry(
+                ifc_file,
+                parent_boundary,
+                exterior_boundary_polygon,
+                face_matrix,
+                unit_scale,
+            )
+            _set_boundary_name(parent_boundary)
+            boundaries.append(parent_boundary)
+            covered_by_face[space_ngon_idx].append(exterior_boundary_polygon)
+
+            boundaries.extend(
+                _process_openings(
+                    ifc_file,
+                    openings_to_process,
+                    face_matrix,
+                    boundary_class,
+                    parent_boundary,
+                    space,
+                    unit_scale,
+                    processed_fillings,
+                    covered_by_face[space_ngon_idx],
+                )
+            )
+
+    boundaries.extend(
+        _fill_face_gaps(
+            ifc_file,
+            space,
+            boundary_class,
+            unit_scale,
+            space_face_polygons,
+            face_matrices,
+            face_matrix_invs,
+            space_face_normals_world,
+            space_verts_local,
+            space_ngons,
+            space_matrix,
+            space_matrix_inv,
+            shapes,
+            element_ngons,
+            covered_by_face,
+            building_elements,
+            all_filling_ids,
+            matched_walls_and_columns,
+        )
+    )
+
+    return boundaries
+
+
+def _match_element_to_space_face(
+    element,
+    shapes,
+    element_ngons,
+    space_matrix,
+    space_matrix_inv,
+    space_verts_l,
+    space_face_polygon,
+    face_matrix_inv,
+    space_face_normal_world,
+):
+    """Match a building element's faces against a single space face.
+
+    :return: A tuple ``(dist_min, matching_polygons, matched_elem_normal)`` with
+        the minimum face distance, the matching boundary polygons and the matched
+        face normal in world space, or ``None`` when the element does not bound
+        this space face.
+    """
+    element_shape = shapes[element.id()]
+    element_matrix = element_shape["matrix"]
+    element_matrix_3x3 = element_matrix[:3, :3]
+    element_matrix_inv = np.linalg.inv(element_matrix)
+
+    element_centroid_world = sb.np_apply_matrix(np.mean(element_shape["verts"], axis=0)[np.newaxis], element_matrix)[0]
+
+    space_centroid = np.mean(space_verts_l, axis=0)
+
+    matching_polygons = []
+    matched_elem_normal = None
+    dist_min = None
+    plane_offset_min = None
+
+    for ngon in element_ngons[element.id()]:
+        elem_verts_l = element_shape["verts"][ngon]
+        elem_face_normal_local = _face_normal(elem_verts_l)
+        if elem_face_normal_local is None:
+            continue
+        elem_face_normal_local = _ensure_outward(
+            elem_face_normal_local, elem_verts_l, element_centroid_world, element_matrix
+        )
+        elem_face_normal_world = element_matrix_3x3 @ elem_face_normal_local
+
+        angle = degrees(acos(max(min(float(np.dot(space_face_normal_world, elem_face_normal_world)), 1), -1)))
+        is_horizontal_face = abs(space_face_normal_world[2]) > 0.5 and abs(elem_face_normal_world[2]) > 0.5
+        is_valid_element = (
+            element.is_a("IfcVirtualElement")
+            or element.is_a("IfcSlab")
+            or element.is_a("IfcWindow")
+            or element.is_a("IfcDoor")
+        )
+        is_anti_parallel = _is_x(angle, 180, tolerance=2)
+        is_parallel = _is_x(angle, 0, tolerance=2)
+
+        if not (is_anti_parallel or (is_horizontal_face and is_parallel and is_valid_element)):
+            continue
+
+        sv_in_elem = sb.np_apply_matrix(space_centroid[np.newaxis], element_matrix_inv @ space_matrix)[0]
+        dist = float(np.dot(sv_in_elem - elem_verts_l[0], elem_face_normal_local))
+        dist_tol = 0.05 if is_horizontal_face else 0.5
+        if abs(dist) > dist_tol:
+            continue
+
+        elem_verts_in_space = sb.np_apply_matrix(elem_verts_l, space_matrix_inv @ element_matrix)
+        face_polygon = _verts_to_polygon(elem_verts_in_space, face_matrix_inv, snap=1e-6)
+        if not face_polygon.is_valid:
+            face_polygon = face_polygon.buffer(0)
+
+        try:
+            gross_boundary_polygon = space_face_polygon.intersection(face_polygon)
+        except shapely.errors.GEOSException:
+            logger.warning(
+                "Skipping invalid geometry for %s (shapely topology error).",
+                element.Name or element.is_a(),
+                exc_info=True,
+            )
+            continue
+
+        if gross_boundary_polygon.is_empty or gross_boundary_polygon.area < 1e-4:
+            continue
+
+        if type(gross_boundary_polygon) == shapely.GeometryCollection:
+            for geom in gross_boundary_polygon.geoms:
+                if type(geom) == shapely.Polygon:
+                    gross_boundary_polygon = geom
+                    break
+
+        if not (isinstance(gross_boundary_polygon, shapely.Polygon) and gross_boundary_polygon.is_valid):
+            continue
+        if gross_boundary_polygon.is_empty:
+            continue
+
+        matching_polygons.append(gross_boundary_polygon)
+        matched_elem_normal = elem_face_normal_world
+        dist_min = abs(dist) if dist_min is None else min(dist_min, abs(dist))
+        space_face_normal = _face_normal(space_verts_l)
+        if space_face_normal is not None:
+            plane_offset = abs(float(np.dot(space_face_normal, elem_verts_in_space[0] - space_verts_l[0])))
+            plane_offset_min = plane_offset if plane_offset_min is None else min(plane_offset_min, plane_offset)
+
+    if not matching_polygons:
+        return None
+    return dist_min, plane_offset_min, matching_polygons, matched_elem_normal
+
+
+def _union_openings_into_parent(exterior_boundary_polygon, openings_to_process):
+    """Union the opening polygons back into the parent boundary polygon.
+
+    Authoring tools may bake openings into the building element mesh, so the
+    parent boundary polygon can be notched where the opening is. Since an inner
+    boundary is supposed to overlap its parent boundary, the openings are
+    unioned back into the parent while the filling gets its own boundary.
+    """
+    for _, _, opening_polygon in openings_to_process:
+        unionised_object = exterior_boundary_polygon.union(opening_polygon)
+        if isinstance(unionised_object, shapely.Polygon):
+            exterior_boundary_polygon = unionised_object
+    return exterior_boundary_polygon
+
+
+def _compute_opening_polygon(ifc_file, opening, face_normal_world, space_matrix_inv, face_matrix_inv):
+    """Project an opening onto the building element face in space-local coordinates.
+
+    :param opening: The IfcOpeningElement to project.
+    :param face_normal_world: The building element face normal in world space.
+    :param space_matrix_inv: Inverse of the space placement matrix.
+    :param face_matrix_inv: The inverse face matrix (for 2D projection).
+    :return: A 2D shapely polygon in space-local coordinates, or None.
+    """
+    settings = ifcopenshell.geom.settings()
+    try:
+        shape = ifcopenshell.geom.create_shape(settings, opening)
+    except Exception:
+        return None
+    opening_verts_l = ifcopenshell.util.shape.get_vertices(shape.geometry)
+    opening_faces = ifcopenshell.util.shape.get_faces(shape.geometry)
+    opening_edges = ifcopenshell.util.shape.get_edges(shape.geometry)
+    opening_matrix = ifcopenshell.util.shape.get_shape_matrix(shape)
+    opening_matrix_3x3 = opening_matrix[:3, :3]
+
+    opening_ngons = ifcopenshell.util.shape.dissolve_faces(
+        opening_verts_l, opening_faces, opening_edges, merge_coplanar=True
+    )
+
+    opening_polygons = []
+    for ngon in opening_ngons:
+        o_verts_l = opening_verts_l[ngon]
+        o_normal_local = _face_normal(o_verts_l)
+        if o_normal_local is None:
+            continue
+        o_normal_world = opening_matrix_3x3 @ o_normal_local
+        angle = degrees(acos(max(min(float(np.dot(o_normal_world, face_normal_world)), 1), -1)))
+        if not _is_x(angle, 180, tolerance=2):
+            continue
+        o_verts_in_space = sb.np_apply_matrix(o_verts_l, space_matrix_inv @ opening_matrix)
+        polygon = _verts_to_polygon(o_verts_in_space, face_matrix_inv)
+        opening_polygons.append(polygon)
+
+    if not opening_polygons:
+        return None
+
+    return shapely.ops.unary_union(opening_polygons)
+
+
+def _process_openings(
+    ifc_file,
+    openings_to_process,
+    face_matrix,
+    boundary_class,
+    parent_boundary,
+    space,
+    unit_scale,
+    processed_fillings: set[int],
+    covered_polygons: list,
+):
+    """Create boundaries for the fillings of openings in a building element.
+
+    :param openings_to_process: Tuples of (opening, filling, opening polygon).
+    :param face_matrix: The face matrix in space-local coordinates (for connection geometry).
+    :param processed_fillings: Set of element IDs that already have opening boundaries.
+    :param covered_polygons: Accumulated boundary polygons used for water-tightness checks.
+    """
+    boundaries = []
+
+    for opening, filling, opening_polygon in openings_to_process:
+        filling_id = filling.id()
+        if filling_id in processed_fillings:
+            continue
+
+        boundary = ifcopenshell.api.root.create_entity(ifc_file, ifc_class=boundary_class)
+        boundary.RelatingSpace = space
+        boundary.RelatedBuildingElement = filling or opening
+
+        # Use the same space-local face_matrix for connection geometry
+        _assign_connection_geometry(
+            ifc_file,
+            boundary,
+            opening_polygon,
+            face_matrix,
+            unit_scale,
+        )
+        if filling:
+            boundary.PhysicalOrVirtualBoundary = "PHYSICAL"
+        else:
+            boundary.PhysicalOrVirtualBoundary = "VIRTUAL"
+        boundary.InternalOrExternalBoundary = parent_boundary.InternalOrExternalBoundary
+        if boundary.is_a() != "IfcRelSpaceBoundary":
+            boundary.ParentBoundary = parent_boundary
+        _set_boundary_name(boundary)
+        processed_fillings.add(filling_id)
+        covered_polygons.append(opening_polygon)
+        boundaries.append(boundary)
+
+    return boundaries
+
+
+def _fill_face_gaps(
+    ifc_file,
+    space,
+    boundary_class,
+    unit_scale,
+    space_face_polygons,
+    face_matrices,
+    face_matrix_invs,
+    space_face_normals_world,
+    space_verts_local,
+    space_ngons,
+    space_matrix,
+    space_matrix_inv,
+    shapes,
+    element_ngons,
+    covered_by_face,
+    building_elements,
+    all_filling_ids,
+    matched_walls_and_columns,
+):
+    """Create boundaries for uncovered parts of space faces to keep them water tight."""
+    boundaries = []
+    for face_idx, space_face_polygon in space_face_polygons.items():
+        covered_polygons = covered_by_face.get(face_idx, [])
+        if not covered_polygons:
+            continue
+        uncovered = space_face_polygon.difference(shapely.ops.unary_union(covered_polygons))
+        if uncovered.is_empty:
+            continue
+        if isinstance(uncovered, shapely.Polygon):
+            fragments = [uncovered]
+        elif isinstance(uncovered, shapely.MultiPolygon):
+            fragments = list(uncovered.geoms)
+        else:
+            continue
+        for fragment in fragments:
+            if fragment.area < 1e-2:
+                continue
+            element = _best_element_for_gap(
+                fragment,
+                space_verts_local[space_ngons[face_idx]],
+                space_matrix,
+                space_matrix_inv,
+                face_matrices[face_idx],
+                face_matrix_invs[face_idx],
+                space_face_normals_world[face_idx],
+                shapes,
+                element_ngons,
+                building_elements,
+                all_filling_ids,
+                matched_walls_and_columns,
+            )
+            if element is None:
+                logger.warning(
+                    "No element found to fill a gap on a face of space %s.",
+                    space.Name or space.is_a(),
+                )
+                continue
+            parent_boundary = ifcopenshell.api.root.create_entity(ifc_file, ifc_class=boundary_class)
+            parent_boundary.PhysicalOrVirtualBoundary = "PHYSICAL"
+            parent_boundary.InternalOrExternalBoundary = "NOTDEFINED"
+            _set_internal_external(parent_boundary, element)
+            parent_boundary.RelatingSpace = space
+            parent_boundary.RelatedBuildingElement = element
+            _assign_connection_geometry(
+                ifc_file,
+                parent_boundary,
+                fragment,
+                face_matrices[face_idx],
+                unit_scale,
+            )
+            _set_boundary_name(parent_boundary)
+            boundaries.append(parent_boundary)
+    return boundaries
+
+
+def _best_element_for_gap(
+    fragment,
+    space_face_verts,
+    space_matrix,
+    space_matrix_inv,
+    face_matrix,
+    face_matrix_inv,
+    space_face_normal_world,
+    shapes,
+    element_ngons,
+    building_elements,
+    all_filling_ids,
+    matched_walls_and_columns,
+):
+    """Find the element most appropriate to cover an uncovered part of a space face."""
+    best = None
+    best_overlap = 0.0
+    space_centroid = np.mean(space_face_verts, axis=0)
+
+    for element in building_elements:
+        if element.id() in all_filling_ids:
+            continue
+        # Walls and columns already bounding this space keep a single boundary;
+        # a gap is therefore filled by a neighbouring element instead.
+        if element.is_a() in ("IfcWall", "IfcColumn") and element.id() in matched_walls_and_columns:
+            continue
+        es = shapes[element.id()]
+        e_matrix = es["matrix"]
+        e_matrix_3x3 = e_matrix[:3, :3]
+        e_matrix_inv = np.linalg.inv(e_matrix)
+        e_centroid_world = sb.np_apply_matrix(np.mean(es["verts"], axis=0)[np.newaxis], e_matrix)[0]
+
+        for ngon in element_ngons[element.id()]:
+            elem_verts_l = es["verts"][ngon]
+            normal_local = _face_normal(elem_verts_l)
+            if normal_local is None:
+                continue
+            normal_local = _ensure_outward(normal_local, elem_verts_l, e_centroid_world, e_matrix)
+            normal_world = e_matrix_3x3 @ normal_local
+            angle = degrees(acos(max(min(float(np.dot(space_face_normal_world, normal_world)), 1), -1)))
+            is_horizontal_face = abs(space_face_normal_world[2]) > 0.5 and abs(normal_world[2]) > 0.5
+            is_valid_element = (
+                element.is_a("IfcVirtualElement")
+                or element.is_a("IfcSlab")
+                or element.is_a("IfcWindow")
+                or element.is_a("IfcDoor")
+            )
+            if not (
+                _is_x(angle, 180, tolerance=2)
+                or (is_horizontal_face and _is_x(angle, 0, tolerance=2) and is_valid_element)
+            ):
+                continue
+            sv_in_elem = sb.np_apply_matrix(space_centroid[np.newaxis], e_matrix_inv @ space_matrix)[0]
+            dist = float(np.dot(sv_in_elem - elem_verts_l[0], normal_local))
+            dist_tol = 0.05 if is_horizontal_face else 0.5
+            if abs(dist) > dist_tol:
+                continue
+            elem_verts_in_space = sb.np_apply_matrix(elem_verts_l, space_matrix_inv @ e_matrix)
+            face_polygon = _verts_to_polygon(elem_verts_in_space, face_matrix_inv, snap=1e-6)
+            if not face_polygon.is_valid:
+                face_polygon = face_polygon.buffer(0)
+            overlap = fragment.intersection(face_polygon).area
+            if overlap > best_overlap:
+                best_overlap = overlap
+                best = element
+
+    if best is not None:
+        return best
+
+    # For gaps at corners between non-parallel faces, fall back to the element
+    # whose plan footprint covers the gap centroid.
+    frag_2d = np.array([[c[0], c[1], 0.0] for c in fragment.exterior.coords])
+    frag_local = sb.np_apply_matrix(frag_2d, face_matrix)
+    frag_world = sb.np_apply_matrix(frag_local, space_matrix)
+    frag_centroid = frag_world.mean(axis=0)
+    is_horizontal_face = abs(space_face_normal_world[2]) > 0.5
+    best = None
+    best_dist = np.inf
+    for element in building_elements:
+        if element.id() in all_filling_ids:
+            continue
+        if is_horizontal_face:
+            if not (element.is_a("IfcSlab") or element.is_a("IfcRoof") or element.is_a("IfcVirtualElement")):
+                continue
+        elif not (element.is_a("IfcWall") or element.is_a("IfcColumn") or element.is_a("IfcVirtualElement")):
+            continue
+        es = shapes[element.id()]
+        world = sb.np_apply_matrix(es["verts"], es["matrix"])
+        elem_xy = shapely.box(world[:, 0].min(), world[:, 1].min(), world[:, 0].max(), world[:, 1].max())
+        if not elem_xy.contains(shapely.Point(frag_centroid[:2])):
+            continue
+        elem_centroid = world.mean(axis=0)
+        dist = float(np.linalg.norm(elem_centroid - frag_centroid))
+        if dist < best_dist:
+            best_dist = dist
+            best = element
+    return best
+
+
+def _face_normal(verts: np.ndarray) -> Optional[np.ndarray]:
+    """Compute the normal of a polygon from its vertices."""
+    if len(verts) < 3:
+        return None
+    for i in range(len(verts) - 2):
+        v0, v1, v2 = verts[i], verts[i + 1], verts[i + 2]
+        cross = np.cross(v1 - v0, v2 - v0)
+        norm = np.linalg.norm(cross)
+        if norm > 1e-8:
+            return cross / norm
+    return None
+
+
+def _face_matrix_from_verts(verts3: np.ndarray) -> np.ndarray:
+    """Build a 4x4 face-local coordinate matrix from 3 vertices."""
+    p1, p2, p3 = verts3[0], verts3[1], verts3[2]
+    z = sb.np_normal([p1, p2, p3])
+    x = sb.np_normalized(p2 - p1)
+    return ifcopenshell.util.placement.a2p(o=p1, z=z, x=x)
+
+
+def _verts_to_polygon(verts: np.ndarray, face_matrix_inv: np.ndarray, snap: float = 0) -> shapely.Polygon:
+    """Project 3D vertices onto a 2D plane and create a shapely Polygon."""
+    verts_2d = sb.np_apply_matrix(verts, face_matrix_inv)[:, :2]
+    if snap:
+        verts_2d = np.round(verts_2d / snap) * snap
+    return shapely.Polygon([tuple(v) for v in verts_2d])
+
+
+def _assign_connection_geometry(
+    ifc_file: ifcopenshell.file,
+    boundary: ifcopenshell.entity_instance,
+    polygon: shapely.Polygon,
+    face_matrix: np.ndarray,
+    unit_scale: float,
+) -> None:
+    """Assign connection geometry to a boundary using the existing API."""
+    location = face_matrix[:3, 3]
+    axis = face_matrix[:3, 2]
+    ref_direction = face_matrix[:3, 0]
+
+    outer_boundary = [list(coord) for coord in polygon.exterior.coords[:-1]]
+    inner_boundaries = [list(interior.coords[:-1]) for interior in polygon.interiors]
+
+    ifcopenshell.api.boundary.assign_connection_geometry(
+        ifc_file,
+        rel_space_boundary=boundary,
+        outer_boundary=outer_boundary,
+        location=location.tolist(),
+        axis=axis.tolist(),
+        ref_direction=ref_direction.tolist(),
+        inner_boundaries=inner_boundaries if inner_boundaries else None,
+        unit_scale=unit_scale,
+    )
+
+
+def _set_internal_external(
+    boundary: ifcopenshell.entity_instance, building_element: ifcopenshell.entity_instance
+) -> None:
+    """Set InternalOrExternalBoundary based on element type and psets."""
+    if building_element.is_a("IfcWall"):
+        is_external = ifcopenshell.util.element.get_pset(building_element, "Pset_WallCommon", "IsExternal")
+        if is_external is True:
+            boundary.InternalOrExternalBoundary = "EXTERNAL"
+        elif is_external is False:
+            boundary.InternalOrExternalBoundary = "INTERNAL"
+    elif building_element.is_a("IfcSlab"):
+        predefined_type = ifcopenshell.util.element.get_predefined_type(building_element)
+        if predefined_type == "BASESLAB":
+            boundary.InternalOrExternalBoundary = "EXTERNAL_EARTH"
+        else:
+            is_external = ifcopenshell.util.element.get_pset(building_element, "Pset_SlabCommon", "IsExternal")
+            if is_external is True:
+                boundary.InternalOrExternalBoundary = "EXTERNAL"
+            elif is_external is False:
+                boundary.InternalOrExternalBoundary = "INTERNAL"
+
+
+def _set_boundary_name(boundary: ifcopenshell.entity_instance) -> None:
+    """Set Name/Description per IFC4x3 convention."""
+    if boundary.is_a("IfcRelSpaceBoundary2ndLevel"):
+        boundary.Name = "2ndLevel"
+        if boundary.CorrespondingBoundary:
+            boundary.Description = "2a"
+        else:
+            boundary.Description = "2b"
+    elif boundary.is_a("IfcRelSpaceBoundary1stLevel"):
+        boundary.Name = "1stLevel"
+
+
+def _ensure_outward(
+    normal_local: np.ndarray,
+    face_verts_l: np.ndarray,
+    entity_centroid_world: np.ndarray,
+    entity_matrix: np.ndarray,
+) -> np.ndarray:
+    """Flip face normal to point away from the entity centroid."""
+    face_centroid_world = sb.np_apply_matrix(np.mean(face_verts_l, axis=0)[np.newaxis], entity_matrix)[0]
+    normal_world = entity_matrix[:3, :3] @ normal_local
+    if np.dot(face_centroid_world - entity_centroid_world, normal_world) < 0:
+        return -normal_local
+    return normal_local
+
+
+def _is_x(value: float, x: float, tolerance: float = 1e-5) -> bool:
+    """Check whether value is within tolerance of x."""
+    return (x + tolerance) > value > (x - tolerance)

src/ifcopenshell-python/ifcopenshell/util/element.py

@@ -2007,3 +2007,24 @@ def get_material_profiles(element: ifcopenshell.entity_instance) -> list[Priorit
         )
         for material_profile in material.MaterialProfiles
     ]
+
+
+def iter_top_connections(
+    element: ifcopenshell.entity_instance,
+) -> Generator[tuple[ifcopenshell.entity_instance, ifcopenshell.entity_instance], None, None]:
+    """Yield ``(connected_element, rel)`` tuples for every
+    ``IfcRelConnectsElements`` with ``Description == "TOP"`` connecting
+    to this element.
+
+    Walks ``element.ConnectedFrom`` because the connecting element (e.g. a
+    slab) is the relating side of the TOP relationship.
+
+    :param element: The IFC element (typically a wall).
+    :return: Generator of ``(connected_element, rel)`` tuples.
+    """
+    for rel in getattr(element, "ConnectedFrom", []) or ():
+        if not rel.is_a("IfcRelConnectsElements") or rel.Description != "TOP":
+            continue
+        connected = rel.RelatingElement
+        if connected is not None:
+            yield connected, rel

src/ifcopenshell-python/ifcopenshell/util/shape.py

@@ -752,3 +752,279 @@ def get_total_edge_length(geometry: W.triangulation) -> float:
     vertices = get_vertices(geometry)
     vertices = vertices[get_edges(geometry)]
     return np.linalg.norm(vertices[:, 1] - vertices[:, 0], axis=1).sum().item()
+
+
+def _extend_line(start: np.ndarray, end: np.ndarray, distance: float) -> tuple[np.ndarray, np.ndarray]:
+    """Extend a line segment by a fixed distance on both ends.
+
+    :param start: (x, y) or (x, y, z) array.
+    :param end: (x, y) or (x, y, z) array.
+    :param distance: Distance to extend on each end.
+    :return: (new_start, new_end) arrays.
+    """
+    direction = end - start
+    norm = np.linalg.norm(direction)
+    if norm == 0:
+        return start, end
+    offset = distance * (direction / norm)
+    return start - offset, end + offset
+
+
+def bisect_mesh_plane_vf(
+    verts: npt.NDArray[np.float64],
+    faces: npt.NDArray[np.int32],
+    plane_z: float,
+    *,
+    precision: int = 3,
+    extend: float = 0.0,
+) -> list:
+    """Intersect a triangulated mesh with a horizontal Z plane.
+
+    All faces are processed at once via numpy broadcasting for performance.
+
+    :param verts: (n, 3) array of vertices in world coordinates.
+    :param faces: (m, 3) array of triangle vertex indices.
+    :param plane_z: Z elevation of the horizontal cutting plane.
+    :param precision: Decimal places to round intersection point coordinates to.
+    :param extend: Distance to extend each segment on both ends, to ensure
+        overlap with neighbouring segments for polygon closure.
+    :return: List of (start_xy, end_xy) tuples where each coordinate is (x, y).
+    """
+    if len(faces) == 0:
+        return []
+    v0 = verts[faces[:, 0]]
+    v1 = verts[faces[:, 1]]
+    v2 = verts[faces[:, 2]]
+    d0 = v0[:, 2] - plane_z
+    d1 = v1[:, 2] - plane_z
+    d2 = v2[:, 2] - plane_z
+    straddle = ~((np.minimum(np.minimum(d0, d1), d2) > 0) | (np.maximum(np.maximum(d0, d1), d2) < 0))
+    if not np.any(straddle):
+        return []
+    idx = np.where(straddle)[0]
+    d0s, d1s, d2s = d0[idx], d1[idx], d2[idx]
+    v0s, v1s, v2s = v0[idx], v1[idx], v2[idx]
+
+    def _edge_intersections(va, vb, da, db):
+        mask = da * db < 0
+        diff = da - db
+        diff = np.where(diff == 0, 1.0, diff)
+        t = np.where(mask, da / diff, 0.0)
+        pts = va + t[:, np.newaxis] * (vb - va)
+        return pts, mask
+
+    p01, m01 = _edge_intersections(v0s, v1s, d0s, d1s)
+    p12, m12 = _edge_intersections(v1s, v2s, d1s, d2s)
+    p20, m20 = _edge_intersections(v2s, v0s, d2s, d0s)
+
+    segments = []
+    for i in range(len(idx)):
+        pts_xy = []
+        for pt, mask in ((p01[i], m01[i]), (p12[i], m12[i]), (p20[i], m20[i])):
+            if mask:
+                pts_xy.append((round(float(pt[0]), precision), round(float(pt[1]), precision)))
+        if len(pts_xy) == 2 and pts_xy[0] != pts_xy[1]:
+            if extend > 0:
+                s, e = _extend_line(np.array(pts_xy[0]), np.array(pts_xy[1]), extend)
+                segments.append((s.tolist(), e.tolist()))
+            else:
+                segments.append(pts_xy)
+    return segments
+
+
+def dissolve_faces(
+    verts: npt.NDArray[np.float64],
+    faces: npt.NDArray[np.int32],
+    edges: npt.NDArray[np.int32],
+    merge_coplanar: bool = False,
+    angle_tolerance: float = 0.017453292519943295,
+) -> list[list[int]]:
+    """Reconstruct polygonal faces from triangulated mesh data.
+
+    Uses the original (pre-triangulation) edges from ``get_edges`` to
+    identify which triangle edges are internal (to be merged) vs external
+    (ngon boundaries). Triangles connected by internal edges are grouped
+    into polygonal faces.
+
+    When ``merge_coplanar`` is True, a second pass merges adjacent ngons
+    whose face normals are parallel within ``angle_tolerance`` radians.
+    This mirrors ``bmesh.ops.dissolve_limit`` behavior where coplanar
+    faces sharing an edge are merged regardless of the original face
+    structure. This is needed when the IFC representation splits a single
+    planar face into multiple faces (e.g. an L-shaped top face split into
+    triangles + quads).
+
+    :param verts: (n, 3) array of vertices.
+    :param faces: (m, 3) array of triangle vertex indices.
+    :param edges: (e, 2) array of original (pre-triangulation) edge vertex
+        indices, as returned by :func:`get_edges`.
+    :param merge_coplanar: If True, merge adjacent coplanar ngons.
+    :param angle_tolerance: Angle in radians for coplanar merge (default 1°).
+    :return: List of polygonal faces, each as an ordered list of vertex indices
+        forming a closed polygon (last vertex connects back to first).
+    """
+    if len(faces) == 0:
+        return []
+    if len(edges) == 0:
+        return [list(f) for f in faces]
+
+    original_edges = {frozenset((int(e[0]), int(e[1]))) for e in edges}
+
+    tri_edges = []
+    for f in faces:
+        tri_edges.append(
+            (
+                frozenset((int(f[0]), int(f[1]))),
+                frozenset((int(f[1]), int(f[2]))),
+                frozenset((int(f[2]), int(f[0]))),
+            )
+        )
+
+    internal_edge_to_tris: dict[frozenset, list[int]] = {}
+    for tri_idx, edges_3 in enumerate(tri_edges):
+        for e in edges_3:
+            if e not in original_edges:
+                internal_edge_to_tris.setdefault(e, []).append(tri_idx)
+
+    parent = list(range(len(faces)))
+
+    def find(x):
+        while parent[x] != x:
+            parent[x] = parent[parent[x]]
+            x = parent[x]
+        return x
+
+    def union(x, y):
+        px, py = find(x), find(y)
+        if px != py:
+            parent[px] = py
+
+    for tri_indices in internal_edge_to_tris.values():
+        if len(tri_indices) == 2:
+            union(tri_indices[0], tri_indices[1])
+
+    ngons: dict[int, list[int]] = {}
+    for tri_idx in range(len(faces)):
+        root = find(tri_idx)
+        ngons.setdefault(root, []).append(tri_idx)
+
+    if merge_coplanar:
+        _merge_coplanar_ngons(ngons, faces, verts, tri_edges, parent, find, union, angle_tolerance)
+
+    result = []
+    for tri_indices in ngons.values():
+        tri_edge_set = set()
+        edge_count: dict[frozenset, int] = {}
+        for tri_idx in tri_indices:
+            for e in tri_edges[tri_idx]:
+                tri_edge_set.add(e)
+                edge_count[e] = edge_count.get(e, 0) + 1
+
+        if merge_coplanar:
+            boundary_edges = [e for e in tri_edge_set if edge_count.get(e, 0) == 1]
+        else:
+            boundary_edges = [e for e in tri_edge_set if e in original_edges]
+
+        if not boundary_edges:
+            result.append(list(faces[tri_indices[0]]))
+            continue
+
+        edge_adjacency: dict[int, int] = {}
+        for e in boundary_edges:
+            v_list = list(e)
+            for tri_idx in tri_indices:
+                f = faces[tri_idx]
+                f_edges = [(int(f[0]), int(f[1])), (int(f[1]), int(f[2])), (int(f[2]), int(f[0]))]
+                for fe in f_edges:
+                    if frozenset(fe) == e:
+                        edge_adjacency[fe[0]] = fe[1]
+                        break
+                else:
+                    continue
+                break
+
+        if not edge_adjacency:
+            result.append(list(faces[tri_indices[0]]))
+            continue
+
+        start = next(iter(edge_adjacency))
+        polygon = [start]
+        current = edge_adjacency[start]
+        while current != start and current in edge_adjacency:
+            polygon.append(current)
+            current = edge_adjacency[current]
+
+        if len(polygon) >= 3:
+            result.append(polygon)
+        else:
+            result.append(list(faces[tri_indices[0]]))
+
+    return result
+
+
+def _merge_coplanar_ngons(
+    ngons: dict[int, list[int]],
+    faces: npt.NDArray[np.int32],
+    verts: npt.NDArray[np.float64],
+    tri_edges: list,
+    parent: list[int],
+    find,
+    union,
+    angle_tolerance: float,
+) -> None:
+    """Merge adjacent ngons whose face normals are parallel within tolerance.
+
+    Modifies ``ngons`` and ``parent`` in place.
+    """
+    from math import acos
+
+    # Compute normal for each ngon
+    ngon_normals: dict[int, np.ndarray] = {}
+    ngon_edge_to_ngons: dict[frozenset, list[int]] = {}
+    ngon_roots = list(ngons.keys())
+
+    for root in ngon_roots:
+        tri_indices = ngons[root]
+        f0 = faces[tri_indices[0]]
+        v0, v1, v2 = verts[f0[0]], verts[f0[1]], verts[f0[2]]
+        edge1 = v1 - v0
+        edge2 = v2 - v0
+        normal = np.cross(edge1, edge2)
+        norm = np.linalg.norm(normal)
+        if norm > 1e-8:
+            normal = normal / norm
+        ngon_normals[root] = normal
+
+        # Collect all edges of this ngon
+        ngon_edges = set()
+        for tri_idx in tri_indices:
+            for e in tri_edges[tri_idx]:
+                ngon_edges.add(e)
+        for e in ngon_edges:
+            ngon_edge_to_ngons.setdefault(e, []).append(root)
+
+    # Find shared edges between different ngons and check coplanarity
+    for edge, root_list in ngon_edge_to_ngons.items():
+        if len(root_list) != 2:
+            continue
+        root_a, root_b = root_list[0], root_list[1]
+        if root_a == root_b:
+            continue
+        # Check if already merged
+        ra, rb = find(root_a), find(root_b)
+        if ra == rb:
+            continue
+        # Compare normals
+        na, nb = ngon_normals[root_a], ngon_normals[root_b]
+        dot = max(min(float(np.dot(na, nb)), 1.0), -1.0)
+        angle = acos(dot)
+        if angle < angle_tolerance:
+            union(root_a, root_b)
+
+    # Rebuild ngons dict with merged groups
+    new_ngons: dict[int, list[int]] = {}
+    for root in ngon_roots:
+        new_root = find(root)
+        new_ngons.setdefault(new_root, []).extend(ngons[root])
+    ngons.clear()
+    ngons.update(new_ngons)

src/ifcopenshell-python/ifcopenshell/util/space.py

@@ -0,0 +1,246 @@
+# IfcOpenShell - IFC toolkit and geometry engine
+# Copyright (C) 2026 Dion Moult <dion@thinkmoult.com>
+#
+# This file is part of IfcOpenShell.
+#
+# IfcOpenShell is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# IfcOpenShell is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with IfcOpenShell.  If not, see <http://www.gnu.org/licenses/>.
+
+"""Blender-independent utilities for space geometry generation.
+
+These functions operate on IFC geometry data (vertices, faces, element
+relationships) without requiring any Blender objects to be loaded. They are
+used by Bonsai's space generation pipeline but can also be used standalone
+for IFC analysis.
+"""
+
+from __future__ import annotations
+
+from typing import Optional, Union
+
+import ifcopenshell
+import ifcopenshell.util.element
+import ifcopenshell.util.shape
+import shapely
+
+BOUNDING_CLASSES = ("IfcWall", "IfcColumn", "IfcMember", "IfcVirtualElement", "IfcPlate")
+HEIGHT_DETECTION_CLASSES = ("IfcSlab", "IfcRoof")
+
+
+def get_boundary_lines(
+    ifc_file: ifcopenshell.file,
+    shapes: dict,
+    cut_z: float,
+    bounding_classes: tuple = BOUNDING_CLASSES,
+) -> tuple[list[shapely.LineString], list[ifcopenshell.entity_instance]]:
+    """Generate boundary lines by bisecting IFC element geometry with a horizontal plane.
+
+    :param ifc_file: The IFC file.
+    :param shapes: Dict of element shapes keyed by element id, as produced by
+        a geometry cache. Each entry must have ``verts`` (n,3 ndarray),
+        ``faces`` (m,3 ndarray), ``bottom_z`` (float), ``top_z`` (float).
+    :param cut_z: Z elevation of the cutting plane in world coordinates.
+    :param bounding_classes: IFC classes to treat as space-bounding elements.
+    :return: ``(boundary_lines, bounding_elements)`` where boundary_lines is a
+        list of shapely LineString segments and bounding_elements is a list of
+        IFC entity instances that intersect the cutting plane.
+    """
+    boundary_lines: list[shapely.LineString] = []
+    bounding_elements: list[ifcopenshell.entity_instance] = []
+
+    for element_id, shape_data in shapes.items():
+        element = ifc_file.by_id(element_id)
+        if not any(element.is_a(cls) for cls in bounding_classes):
+            continue
+        if cut_z <= shape_data["bottom_z"] or cut_z >= shape_data["top_z"]:
+            continue
+        bounding_elements.append(element)
+        segments = ifcopenshell.util.shape.bisect_mesh_plane_vf(
+            shape_data["verts"], shape_data["faces"], cut_z, precision=3, extend=0.05
+        )
+        for start, end in segments:
+            boundary_lines.append(shapely.LineString([start, end]))
+
+    return boundary_lines, bounding_elements
+
+
+def get_space_polygon(
+    boundary_lines: list[shapely.LineString],
+    x: float,
+    y: float,
+) -> tuple[Union[shapely.Polygon, str], list]:
+    """Assemble boundary lines into closed polygons and find the one containing (x, y).
+
+    :param boundary_lines: List of shapely LineString segments forming a planar graph.
+    :param x: X coordinate of the point to test.
+    :param y: Y coordinate of the point to test.
+    :return: ``(polygon, [])`` on success, or ``("NO POLYGONS FOUND", [])`` /
+        ``("NO POLYGON FOR POINT", [])`` on failure. The second element is
+        reserved for bounding elements (returned by the caller from
+        :func:`get_boundary_lines`).
+    """
+    unioned = shapely.union_all(shapely.GeometryCollection(boundary_lines))
+    closed_polygons = shapely.polygonize(unioned.geoms)
+    if not closed_polygons:
+        return "NO POLYGONS FOUND", []
+    for polygon in closed_polygons.geoms:
+        if shapely.contains_xy(polygon, x, y):
+            return shapely.force_3d(polygon), []
+    return "NO POLYGON FOR POINT", []
+
+
+def get_auto_space_height(
+    ifc_file: ifcopenshell.file,
+    shapes: dict,
+    space_polygon: shapely.Polygon,
+    base_z: float,
+    bounding_walls: list[ifcopenshell.entity_instance],
+) -> Optional[float]:
+    """Auto-detect space height from elements above using IFC geometry.
+
+    Detection priority:
+    1. ``IfcRelConnectsElements`` (TOP) connections on bounding walls
+    2. ``IfcSlab`` / ``IfcRoof`` elements above with XY overlap to the space polygon
+    3. Minimum wall top Z of bounding walls
+
+    :param ifc_file: The IFC file.
+    :param shapes: Dict of element shapes keyed by element id (see :func:`get_boundary_lines`).
+    :param space_polygon: The space footprint polygon in world XY.
+    :param base_z: The space's base Z in world coordinates.
+    :param bounding_walls: List of IFC wall elements bounding the space.
+    :return: Detected height in meters, or ``None`` if nothing found.
+    """
+    height = get_height_from_top_connections(ifc_file, shapes, bounding_walls, base_z, space_polygon)
+    if height is not None and height > 0:
+        return height
+
+    height = get_height_from_elements_above(ifc_file, shapes, space_polygon, base_z)
+    if height is not None and height > 0:
+        return height
+
+    height = get_height_from_wall_tops(shapes, bounding_walls, base_z)
+    if height is not None and height > 0:
+        return height
+
+    return None
+
+
+def get_height_from_top_connections(
+    ifc_file: ifcopenshell.file,
+    shapes: dict,
+    bounding_walls: list[ifcopenshell.entity_instance],
+    base_z: float,
+    space_polygon: shapely.Polygon,
+) -> Optional[float]:
+    """Find the lowest bottom face of elements connected to bounding walls via IfcRelConnectsElements(TOP).
+
+    :param ifc_file: The IFC file.
+    :param shapes: Dict of element shapes keyed by element id.
+    :param bounding_walls: List of IFC wall elements.
+    :param base_z: The space's base Z in world coordinates.
+    :param space_polygon: The space footprint polygon in world XY.
+    :return: Height in meters, or ``None``.
+    """
+    lowest_min_z: Optional[float] = None
+    for wall_element in bounding_walls:
+        for connected_element, _rel in ifcopenshell.util.element.iter_top_connections(wall_element):
+            if not (connected_element.is_a("IfcSlab") or connected_element.is_a("IfcRoof")):
+                continue
+            shape_data = shapes.get(connected_element.id())
+            if not shape_data:
+                continue
+            min_z = shape_data["bottom_z"]
+            if min_z <= base_z:
+                continue
+            verts = shape_data["verts"]
+            element_box = shapely.box(
+                float(verts[:, 0].min()),
+                float(verts[:, 1].min()),
+                float(verts[:, 0].max()),
+                float(verts[:, 1].max()),
+            )
+            if not element_box.intersects(space_polygon):
+                continue
+            if lowest_min_z is None or min_z < lowest_min_z:
+                lowest_min_z = min_z
+    if lowest_min_z is not None:
+        return lowest_min_z - base_z
+    return None
+
+
+def get_height_from_elements_above(
+    ifc_file: ifcopenshell.file,
+    shapes: dict,
+    space_polygon: shapely.Polygon,
+    base_z: float,
+    height_classes: tuple = HEIGHT_DETECTION_CLASSES,
+) -> Optional[float]:
+    """Find the lowest IfcSlab / IfcRoof above whose XY bbox overlaps the space polygon.
+
+    :param ifc_file: The IFC file.
+    :param shapes: Dict of element shapes keyed by element id.
+    :param space_polygon: The space footprint polygon in world XY.
+    :param base_z: The space's base Z in world coordinates.
+    :param height_classes: IFC classes to consider as ceiling elements.
+    :return: Height in meters, or ``None``.
+    """
+    lowest_min_z: Optional[float] = None
+    for ifc_class in height_classes:
+        for element in ifc_file.by_type(ifc_class):
+            shape_data = shapes.get(element.id())
+            if not shape_data:
+                continue
+            min_z = shape_data["bottom_z"]
+            if min_z <= base_z:
+                continue
+            verts = shape_data["verts"]
+            element_box = shapely.box(
+                float(verts[:, 0].min()),
+                float(verts[:, 1].min()),
+                float(verts[:, 0].max()),
+                float(verts[:, 1].max()),
+            )
+            if not element_box.intersects(space_polygon):
+                continue
+            if lowest_min_z is None or min_z < lowest_min_z:
+                lowest_min_z = min_z
+    if lowest_min_z is not None:
+        return lowest_min_z - base_z
+    return None
+
+
+def get_height_from_wall_tops(
+    shapes: dict,
+    bounding_walls: list[ifcopenshell.entity_instance],
+    base_z: float,
+) -> Optional[float]:
+    """Find the minimum wall top Z among bounding walls.
+
+    :param shapes: Dict of element shapes keyed by element id.
+    :param bounding_walls: List of IFC wall elements.
+    :param base_z: The space's base Z in world coordinates.
+    :return: Height in meters, or ``None``.
+    """
+    lowest_top_z: Optional[float] = None
+    for wall_element in bounding_walls:
+        shape_data = shapes.get(wall_element.id())
+        if not shape_data:
+            continue
+        max_z = shape_data["top_z"]
+        if max_z <= base_z:
+            continue
+        if lowest_top_z is None or max_z < lowest_top_z:
+            lowest_top_z = max_z
+    if lowest_top_z is not None:
+        return lowest_top_z - base_z
+    return None

src/ifcopenshell-python/test/util/test_boundary.py

@@ -0,0 +1,403 @@
+# IfcOpenShell - IFC toolkit and geometry engine
+# Copyright (C) 2026 Dion Moult <dion@thinkmoult.com>
+#
+# This file is part of IfcOpenShell.
+#
+# IfcOpenShell is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# IfcOpenShell is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with IfcOpenShell.  If not, see <http://www.gnu.org/licenses/>.
+
+import os
+from collections import Counter
+
+import numpy as np
+import pytest
+import shapely
+
+import ifcopenshell.api.geometry
+import ifcopenshell.api.root
+import ifcopenshell.geom
+import ifcopenshell.util.boundary as subject
+import ifcopenshell.util.shape
+import test.bootstrap
+
+
+def _add_extruded_body(ifc_file, element, coords_2d, depth, z_offset=0.0):
+    """Add a body representation (extruded polyline) to an element."""
+    if not ifc_file.by_type("IfcProject"):
+        ifcopenshell.api.root.create_entity(ifc_file, ifc_class="IfcProject")
+    ctx = ifc_file.createIfcGeometricRepresentationContext(
+        ContextType="Model",
+        CoordinateSpaceDimension=3,
+        Precision=1e-5,
+        WorldCoordinateSystem=ifc_file.createIfcAxis2Placement3D(
+            ifc_file.createIfcCartesianPoint((0.0, 0.0, 0.0)),
+            ifc_file.createIfcDirection((0.0, 0.0, 1.0)),
+            ifc_file.createIfcDirection((1.0, 0.0, 0.0)),
+        ),
+    )
+    sub_ctx = ifc_file.createIfcGeometricRepresentationSubContext(
+        ContextIdentifier="Body",
+        ContextType="Model",
+        ParentContext=ctx,
+        TargetView="MODEL_VIEW",
+    )
+    pts = [ifc_file.createIfcCartesianPoint((float(x), float(y))) for x, y in coords_2d]
+    polyline = ifc_file.createIfcPolyline(pts)
+    profile = ifc_file.create_entity("IfcArbitraryClosedProfileDef", ProfileType="CURVE", OuterCurve=polyline)
+    placement = ifc_file.createIfcAxis2Placement3D(
+        ifc_file.createIfcCartesianPoint((0.0, 0.0, z_offset)),
+        ifc_file.createIfcDirection((0.0, 0.0, 1.0)),
+        ifc_file.createIfcDirection((1.0, 0.0, 0.0)),
+    )
+    direction = ifc_file.createIfcDirection((0.0, 0.0, 1.0))
+    solid = ifc_file.createIfcExtrudedAreaSolid(profile, placement, direction, depth)
+    rep = ifc_file.create_entity(
+        "IfcShapeRepresentation",
+        ContextOfItems=sub_ctx,
+        RepresentationIdentifier="Body",
+        RepresentationType="SweptSolid",
+        Items=[solid],
+    )
+    ifcopenshell.api.geometry.assign_representation(ifc_file, product=element, representation=rep)
+
+
+def _build_shapes_dict(ifc_file, elements):
+    """Build a shapes dict as expected by ifcopenshell.util.boundary."""
+    settings = ifcopenshell.geom.settings()
+    settings.set("disable-opening-subtractions", True)
+    shapes = {}
+    for element in elements:
+        shape = ifcopenshell.geom.create_shape(settings, element)
+        shapes[element.id()] = {
+            "verts": ifcopenshell.util.shape.get_vertices(shape.geometry),
+            "faces": ifcopenshell.util.shape.get_faces(shape.geometry),
+            "edges": ifcopenshell.util.shape.get_edges(shape.geometry),
+            "matrix": ifcopenshell.util.shape.get_shape_matrix(shape),
+        }
+    return shapes
+
+
+def _build_shapes_dict_from_iterator(ifc_file):
+    """Build a shapes dict for all products in a file (excluding openings)."""
+    settings = ifcopenshell.geom.settings()
+    settings.set("disable-opening-subtractions", True)
+    shapes = {}
+    iterator = ifcopenshell.geom.iterator(settings, ifc_file)
+    if iterator.initialize():
+        while True:
+            shape = iterator.get()
+            element = ifc_file.by_id(shape.id)
+            if not element.is_a("IfcOpeningElement"):
+                shapes[shape.id] = {
+                    "verts": ifcopenshell.util.shape.get_vertices(shape.geometry),
+                    "faces": ifcopenshell.util.shape.get_faces(shape.geometry),
+                    "edges": ifcopenshell.util.shape.get_edges(shape.geometry),
+                    "matrix": ifcopenshell.util.shape.get_shape_matrix(shape),
+                }
+            if not iterator.next():
+                break
+    return shapes
+
+
+def _boundaries_for(boundaries, element):
+    return [b for b in boundaries if b.RelatedBuildingElement == element]
+
+
+def _boundary_inner_count(boundary):
+    surface = boundary.ConnectionGeometry.SurfaceOnRelatingElement if boundary.ConnectionGeometry else None
+    if surface and surface.is_a("IfcCurveBoundedPlane") and surface.InnerBoundaries:
+        return len(surface.InnerBoundaries)
+    return 0
+
+
+def _outer_boundary_area(boundary):
+    surface = boundary.ConnectionGeometry.SurfaceOnRelatingElement
+    points = [(p.Coordinates[0], p.Coordinates[1]) for p in surface.OuterBoundary.Points]
+    area = 0.0
+    for (x1, y1), (x2, y2) in zip(points, points[1:]):
+        area += x1 * y2 - x2 * y1
+    return 0.5 * abs(area)
+
+
+def _boundary_polygon_3d(boundary):
+    """The boundary outer boundary as world-space 3D points."""
+    surface = boundary.ConnectionGeometry.SurfaceOnRelatingElement
+    position = surface.BasisSurface.Position
+    origin = np.array(position.Location.Coordinates, dtype=float)
+    z = np.array(position.Axis.DirectionRatios if position.Axis else [0, 0, 1], dtype=float)
+    x = np.array(position.RefDirection.DirectionRatios if position.RefDirection else [1, 0, 0], dtype=float)
+    y = np.cross(z, x)
+    points = np.array([[p.Coordinates[0], p.Coordinates[1]] for p in surface.OuterBoundary.Points])
+    return origin + points[:, 0, None] * x + points[:, 1, None] * y
+
+
+def _boundary_polygon_in_plane(boundary, reference=None):
+    """The boundary polygon projected onto the reference boundary plane."""
+    reference = reference or boundary
+    surface = reference.ConnectionGeometry.SurfaceOnRelatingElement
+    position = surface.BasisSurface.Position
+    origin = np.array(position.Location.Coordinates, dtype=float)
+    z = np.array(position.Axis.DirectionRatios if position.Axis else [0, 0, 1], dtype=float)
+    x = np.array(position.RefDirection.DirectionRatios if position.RefDirection else [1, 0, 0], dtype=float)
+    y = np.cross(z, x)
+    points = _boundary_polygon_3d(boundary) - origin
+    coords = [(float(p @ x), float(p @ y)) for p in points]
+    return shapely.Polygon(coords)
+
+
+def _add_wall_with_window(ifc_file):
+    """Add a space bounded by a wall with a fully interior opening filled by a window."""
+    space = ifcopenshell.api.root.create_entity(ifc_file, ifc_class="IfcSpace")
+    wall = ifcopenshell.api.root.create_entity(ifc_file, ifc_class="IfcWall")
+    opening_element = ifcopenshell.api.root.create_entity(ifc_file, ifc_class="IfcOpeningElement")
+    window = ifcopenshell.api.root.create_entity(ifc_file, ifc_class="IfcWindow")
+    _add_extruded_body(ifc_file, space, [[-5, -5], [5, -5], [5, 5], [-5, 5]], 3.0)
+    _add_extruded_body(ifc_file, wall, [[-5, 5], [5, 5], [5, 5.5], [-5, 5.5]], 3.0)
+    _add_extruded_body(ifc_file, opening_element, [[-2, 5], [2, 5], [2, 5.5], [-2, 5.5]], 1.8, z_offset=0.6)
+    _add_extruded_body(ifc_file, window, [[-2, 4.5], [2, 4.5], [2, 5.5], [-2, 5.5]], 1.5, z_offset=0.75)
+    ifc_file.createIfcRelVoidsElement(RelatingBuildingElement=wall, RelatedOpeningElement=opening_element)
+    ifc_file.createIfcRelFillsElement(RelatingOpeningElement=opening_element, RelatedBuildingElement=window)
+    return space, wall, window
+
+
+def _add_roof_with_skylight(ifc_file):
+    """Add a space with a roof pierced by an opening covered by a skylight window.
+
+    The window is deliberately not related through IfcRelFillsElement to exercise
+    the geometric detection of fillings.
+    """
+    space = ifcopenshell.api.root.create_entity(ifc_file, ifc_class="IfcSpace")
+    roof = ifcopenshell.api.root.create_entity(ifc_file, ifc_class="IfcRoof")
+    opening_element = ifcopenshell.api.root.create_entity(ifc_file, ifc_class="IfcOpeningElement")
+    window = ifcopenshell.api.root.create_entity(ifc_file, ifc_class="IfcWindow")
+    _add_extruded_body(ifc_file, space, [[-5, -5], [5, -5], [5, 5], [-5, 5]], 3.0)
+    _add_extruded_body(ifc_file, roof, [[-5, -5], [5, -5], [5, 5], [-5, 5]], 0.5, z_offset=3.0)
+    _add_extruded_body(ifc_file, opening_element, [[-1, -1], [1, -1], [1, 1], [-1, 1]], 0.8, z_offset=2.8)
+    _add_extruded_body(ifc_file, window, [[-1, -1], [1, -1], [1, 1], [-1, 1]], 0.2, z_offset=3.5)
+    ifc_file.createIfcRelVoidsElement(RelatingBuildingElement=roof, RelatedOpeningElement=opening_element)
+    return space, roof, window
+
+
+def _external_earth_ifczip():
+    return os.path.join(
+        os.path.dirname(__file__),
+        "..",
+        "IfcRelSpaceBoundary_TestFiles",
+        "IfcRelSpaceBoundary2ndLevel",
+        "ExternalEarth_R20_IFC4.ifczip",
+    )
+
+
+def _over_splitted_ifc():
+    return os.path.join(
+        os.path.dirname(__file__),
+        "..",
+        "IfcRelSpaceBoundary_TestFiles",
+        "IfcRelSpaceBoundary2ndLevel",
+        "OverSplitted_R20_IFC2X3.ifc",
+    )
+
+
+def _small_house_ifczip():
+    return os.path.join(
+        os.path.dirname(__file__),
+        "..",
+        "IfcRelSpaceBoundary_TestFiles",
+        "IfcRelSpaceBoundary2ndLevel",
+        "SmallHouse_BB_IFC4.ifczip",
+    )
+
+
+def _triangle_ifczip():
+    return os.path.join(
+        os.path.dirname(__file__),
+        "..",
+        "IfcRelSpaceBoundary_TestFiles",
+        "IfcRelSpaceBoundary2ndLevel",
+        "Triangle_BB_IFC4.ifczip",
+    )
+
+
+def _boundary_element_counts(ifc_file, space_id):
+    space = ifc_file.by_id(space_id)
+    counts = Counter()
+    for boundary in space.BoundedBy or []:
+        if boundary.RelatedBuildingElement:
+            counts[boundary.RelatedBuildingElement.id()] += 1
+    return counts
+
+
+class TestAutoGenerateBoundaries(test.bootstrap.IFC4):
+    def test_no_building_elements_returns_error(self):
+        space = ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcSpace")
+        _add_extruded_body(self.file, space, [[-5, -5], [5, -5], [5, 5], [-5, 5]], 3.0)
+        shapes = _build_shapes_dict(self.file, [space])
+        result = subject.auto_generate_boundaries(self.file, space, shapes, "IfcRelSpaceBoundary")
+        assert isinstance(result, str)
+        assert "No building elements" in result
+
+    def test_space_not_in_shapes_returns_error(self):
+        space = ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcSpace")
+        result = subject.auto_generate_boundaries(self.file, space, {}, "IfcRelSpaceBoundary")
+        assert isinstance(result, str)
+        assert "not found" in result.lower()
+
+    def test_generates_boundary_for_adjacent_wall(self):
+        space = ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcSpace")
+        wall = ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcWall")
+        _add_extruded_body(self.file, space, [[-5, -5], [5, -5], [5, 5], [-5, 5]], 3.0)
+        _add_extruded_body(self.file, wall, [[-5, 5], [5, 5], [5, 5.2], [-5, 5.2]], 3.0)
+        shapes = _build_shapes_dict(self.file, [space, wall])
+        result = subject.auto_generate_boundaries(self.file, space, shapes, "IfcRelSpaceBoundary")
+        assert isinstance(result, list)
+        assert len(result) >= 1
+        boundary = result[0]
+        assert boundary.RelatingSpace == space
+        assert boundary.RelatedBuildingElement == wall
+        assert boundary.PhysicalOrVirtualBoundary == "PHYSICAL"
+
+    def test_wall_boundary_with_window_has_no_inner_boundary(self):
+        space, wall, window = _add_wall_with_window(self.file)
+        shapes = _build_shapes_dict(self.file, [space, wall, window])
+        result = subject.auto_generate_boundaries(self.file, space, shapes, "IfcRelSpaceBoundary2ndLevel")
+        assert isinstance(result, list)
+        wall_boundaries = _boundaries_for(result, wall)
+        assert len(wall_boundaries) == 1
+        assert _boundary_inner_count(wall_boundaries[0]) == 0
+        assert _outer_boundary_area(wall_boundaries[0]) == pytest.approx(30.0, abs=1e-3)
+        window_boundaries = _boundaries_for(result, window)
+        assert len(window_boundaries) == 1
+        assert window_boundaries[0].ParentBoundary == wall_boundaries[0]
+
+    def test_roof_boundary_with_skylight_has_no_inner_boundary(self):
+        space, roof, window = _add_roof_with_skylight(self.file)
+        shapes = _build_shapes_dict(self.file, [space, roof, window])
+        result = subject.auto_generate_boundaries(self.file, space, shapes, "IfcRelSpaceBoundary2ndLevel")
+        assert isinstance(result, list)
+        roof_boundaries = _boundaries_for(result, roof)
+        assert len(roof_boundaries) == 1
+        assert _boundary_inner_count(roof_boundaries[0]) == 0
+        assert _outer_boundary_area(roof_boundaries[0]) == pytest.approx(100.0, abs=1e-3)
+        window_boundaries = _boundaries_for(result, window)
+        assert len(window_boundaries) == 1
+        assert window_boundaries[0].ParentBoundary == roof_boundaries[0]
+
+    def test_wall_boundary_with_unfilled_opening_has_no_inner_boundary(self):
+        space = ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcSpace")
+        wall = ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcWall")
+        opening_element = ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcOpeningElement")
+        _add_extruded_body(self.file, space, [[-5, -5], [5, -5], [5, 5], [-5, 5]], 3.0)
+        _add_extruded_body(self.file, wall, [[-5, 5], [5, 5], [5, 5.5], [-5, 5.5]], 3.0)
+        _add_extruded_body(self.file, opening_element, [[-2, 5], [2, 5], [2, 5.5], [-2, 5.5]], 1.8, z_offset=0.6)
+        self.file.createIfcRelVoidsElement(RelatingBuildingElement=wall, RelatedOpeningElement=opening_element)
+        shapes = _build_shapes_dict(self.file, [space, wall])
+        result = subject.auto_generate_boundaries(self.file, space, shapes, "IfcRelSpaceBoundary2ndLevel")
+        assert isinstance(result, list)
+        wall_boundaries = _boundaries_for(result, wall)
+        assert len(wall_boundaries) == 1
+        assert _boundary_inner_count(wall_boundaries[0]) == 0
+
+    def test_openings_are_unioned_into_parent_boundary(self):
+        # Some authoring tools bake the opening into the building element mesh,
+        # leaving a notch in the gross boundary polygon. The parent boundary must
+        # union the opening back in so it overlaps its own inner boundary.
+        wall_face = shapely.Polygon([(-5, 5), (5, 5), (5, 5.5), (2, 5.5), (2, 5), (-2, 5), (-2, 5.5), (-5, 5.5)])
+        window = shapely.Polygon([(-2, 5), (2, 5), (2, 5.5), (-2, 5.5)])
+        assert wall_face.area == pytest.approx(3.0, abs=1e-9)
+        parent = subject._union_openings_into_parent(wall_face, [("opening", "window", window)])
+        assert isinstance(parent, shapely.Polygon)
+        assert parent.area == pytest.approx(5.0, abs=1e-9)
+        assert parent.contains(window)
+
+    def test_external_earth_boundaries(self):
+        ifczip = _external_earth_ifczip()
+        if not os.path.exists(ifczip):
+            pytest.skip("IfcRelSpaceBoundary_TestFiles submodule is not checked out")
+        ifc_file = ifcopenshell.open(ifczip)
+        shapes = _build_shapes_dict_from_iterator(ifc_file)
+        for space_id, expected_counts in [
+            (182, {996: 1, 1122: 1, 1235: 2, 1288: 1, 2970: 1, 3071: 1, 3140: 1, 3214: 1, 3435: 1, 3605: 1, 3719: 1}),
+            (440, {1122: 1, 3140: 1, 3214: 1, 3493: 1, 3605: 1, 3640: 1, 3669: 1, 3719: 1}),
+            (628, {3140: 1, 3838: 1, 3927: 1, 3980: 2, 4033: 1, 4086: 1, 4139: 1, 4199: 1}),
+        ]:
+            copy = ifcopenshell.file.from_string(ifc_file.wrapped_data.to_string())
+            new_space = copy.by_id(space_id)
+            result = subject.auto_generate_boundaries(
+                copy, new_space, shapes=shapes, boundary_class="IfcRelSpaceBoundary2ndLevel"
+            )
+            assert _boundary_element_counts(copy, space_id) == expected_counts
+            for boundary in result:
+                assert _boundary_inner_count(boundary) == 0
+                if boundary.ParentBoundary:
+                    parent_polygon = _boundary_polygon_in_plane(boundary.ParentBoundary)
+                    child_polygon = _boundary_polygon_in_plane(boundary, reference=boundary.ParentBoundary)
+                    assert child_polygon.intersection(parent_polygon).area == pytest.approx(
+                        child_polygon.area, abs=1e-2
+                    )
+            if space_id in (182, 440):
+                roof_boundaries = _boundaries_for(result, copy.by_id(3214))
+                assert len(roof_boundaries) == 1
+                skylight = [b for b in result if b.RelatedBuildingElement.id() in (3435, 3640)]
+                assert len(skylight) == 1
+                assert skylight[0].ParentBoundary == roof_boundaries[0]
+
+    def test_over_splitted_roof_keeps_shaft_opening(self):
+        # Space 251 of OverSplitted_R20_IFC2X3.ifc has an L-shaped ceiling
+        # pierced by a shaft opening. The ceiling boundary must keep the
+        # opening as an inner boundary, and each shaft wall must get exactly
+        # one boundary instead of fragmented partial + gap boundaries.
+        ifc_path = _over_splitted_ifc()
+        if not os.path.exists(ifc_path):
+            pytest.skip("IfcRelSpaceBoundary_TestFiles submodule is not checked out")
+        ifc_file = ifcopenshell.open(ifc_path)
+        shapes = _build_shapes_dict_from_iterator(ifc_file)
+        copy = ifcopenshell.file.from_string(ifc_file.wrapped_data.to_string())
+        result = subject.auto_generate_boundaries(
+            copy, copy.by_id(251), shapes=shapes, boundary_class="IfcRelSpaceBoundary"
+        )
+        assert isinstance(result, list)
+        assert len(result) == 17
+        roof_boundaries = _boundaries_for(result, copy.by_id(5248))
+        assert len(roof_boundaries) == 1
+        assert _boundary_inner_count(roof_boundaries[0]) == 1
+        surface = roof_boundaries[0].ConnectionGeometry.SurfaceOnRelatingElement
+        outer = [(p.Coordinates[0], p.Coordinates[1]) for p in surface.OuterBoundary.Points]
+        inner = [[(p.Coordinates[0], p.Coordinates[1]) for p in ib.Points] for ib in surface.InnerBoundaries]
+        assert shapely.Polygon(outer, inner).area == pytest.approx(9.962, abs=1e-3)
+        for wall_id in (5832, 5877, 5922, 5967):
+            assert len(_boundaries_for(result, copy.by_id(wall_id))) == 1
+
+    def test_small_house_boundaries(self):
+        ifc_path = _small_house_ifczip()
+        if not os.path.exists(ifc_path):
+            pytest.skip("IfcRelSpaceBoundary_TestFiles submodule is not checked out")
+        ifc_file = ifcopenshell.open(ifc_path)
+        shapes = _build_shapes_dict_from_iterator(ifc_file)
+        for space_id, expected_total in [(1692, 8), (4356, 8), (4380, 13), (6185, 1)]:
+            copy = ifcopenshell.file.from_string(ifc_file.wrapped_data.to_string())
+            result = subject.auto_generate_boundaries(
+                copy, copy.by_id(space_id), shapes=shapes, boundary_class="IfcRelSpaceBoundary2ndLevel"
+            )
+            assert len(result) == expected_total
+
+    def test_triangle_boundaries(self):
+        ifc_path = _triangle_ifczip()
+        if not os.path.exists(ifc_path):
+            pytest.skip("IfcRelSpaceBoundary_TestFiles submodule is not checked out")
+        ifc_file = ifcopenshell.open(ifc_path)
+        shapes = _build_shapes_dict_from_iterator(ifc_file)
+        copy = ifcopenshell.file.from_string(ifc_file.wrapped_data.to_string())
+        result = subject.auto_generate_boundaries(
+            copy, copy.by_id(1573), shapes=shapes, boundary_class="IfcRelSpaceBoundary2ndLevel"
+        )
+        assert len(result) == 6

src/ifcopenshell-python/test/util/test_element.py

@@ -1392,3 +1392,47 @@ class TestCopyDeepIFC4(test.bootstrap.IFC4):
         element2 = subject.copy_deep(self.file, element)
         assert element2.Segments[0][0] == (1, 2)
         assert element2.Segments[1][0] == (3, 4)
+
+
+class TestIterTopConnections(test.bootstrap.IFC4):
+    def test_yields_top_connected_element(self):
+        wall = ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcWall")
+        slab = ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcSlab")
+        rel = self.file.createIfcRelConnectsElements(
+            GlobalId=ifcopenshell.guid.new(),
+            RelatingElement=slab,
+            RelatedElement=wall,
+            Description="TOP",
+        )
+        results = list(subject.iter_top_connections(wall))
+        assert len(results) == 1
+        assert results[0][0] == slab
+        assert results[0][1] == rel
+
+    def test_returns_empty_when_no_connections(self):
+        wall = ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcWall")
+        assert list(subject.iter_top_connections(wall)) == []
+
+    def test_filters_non_top_description(self):
+        wall = ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcWall")
+        slab = ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcSlab")
+        self.file.createIfcRelConnectsElements(
+            GlobalId=ifcopenshell.guid.new(),
+            RelatingElement=slab,
+            RelatedElement=wall,
+            Description="BOTTOM",
+        )
+        assert list(subject.iter_top_connections(wall)) == []
+
+    def test_filters_non_rel_connects_elements(self):
+        wall = ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcWall")
+        slab = ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcSlab")
+        self.file.createIfcRelConnectsPathElements(
+            GlobalId=ifcopenshell.guid.new(),
+            RelatingElement=slab,
+            RelatedElement=wall,
+            Description="ATPATH",
+            RelatingConnectionType="ATPATH",
+            RelatedConnectionType="ATPATH",
+        )
+        assert list(subject.iter_top_connections(wall)) == []

src/ifcopenshell-python/test/util/test_shape.py

@@ -0,0 +1,139 @@
+# IfcOpenShell - IFC toolkit and geometry engine
+# Copyright (C) 2026 Dion Moult <dion@thinkmoult.com>
+#
+# This file is part of IfcOpenShell.
+#
+# IfcOpenShell is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# IfcOpenShell is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with IfcOpenShell.  If not, see <http://www.gnu.org/licenses/>.
+
+import numpy as np
+
+import ifcopenshell.util.shape as subject
+
+
+def _cube_verts_faces(size=2.0, z_offset=0.0):
+    """Build a triangulated cube as (verts, faces) numpy arrays."""
+    s = size / 2
+    verts = np.array(
+        [
+            [-s, -s, -s + z_offset],
+            [s, -s, -s + z_offset],
+            [s, s, -s + z_offset],
+            [-s, s, -s + z_offset],
+            [-s, -s, s + z_offset],
+            [s, -s, s + z_offset],
+            [s, s, s + z_offset],
+            [-s, s, s + z_offset],
+        ],
+        dtype=np.float64,
+    )
+    faces = np.array(
+        [
+            [0, 1, 2],
+            [0, 2, 3],
+            [4, 6, 5],
+            [4, 7, 6],
+            [0, 4, 5],
+            [0, 5, 1],
+            [1, 5, 6],
+            [1, 6, 2],
+            [2, 6, 7],
+            [2, 7, 3],
+            [3, 7, 4],
+            [3, 4, 0],
+        ],
+        dtype=np.int32,
+    )
+    return verts, faces
+
+
+class TestBisectMeshPlaneVf:
+    def test_bisect_at_mid_height(self):
+        verts, faces = _cube_verts_faces(size=2.0)
+        segments = subject.bisect_mesh_plane_vf(verts, faces, plane_z=0.0)
+        assert len(segments) >= 4
+        for start, end in segments:
+            assert len(start) == 2
+            assert len(end) == 2
+
+    def test_bisect_above_mesh_returns_empty(self):
+        verts, faces = _cube_verts_faces(size=2.0)
+        segments = subject.bisect_mesh_plane_vf(verts, faces, plane_z=10.0)
+        assert segments == []
+
+    def test_bisect_below_mesh_returns_empty(self):
+        verts, faces = _cube_verts_faces(size=2.0)
+        segments = subject.bisect_mesh_plane_vf(verts, faces, plane_z=-10.0)
+        assert segments == []
+
+    def test_bisect_with_extend(self):
+        verts, faces = _cube_verts_faces(size=2.0)
+        segments_no_extend = subject.bisect_mesh_plane_vf(verts, faces, plane_z=0.0, extend=0.0)
+        segments_extend = subject.bisect_mesh_plane_vf(verts, faces, plane_z=0.0, extend=0.05)
+        assert len(segments_extend) == len(segments_no_extend)
+        for (s_ext, e_ext), (s_no, e_no) in zip(segments_extend, segments_no_extend):
+            assert abs(s_ext[0] - s_no[0]) >= 0.04 or abs(s_ext[1] - s_no[1]) >= 0.04
+
+    def test_bisect_empty_faces(self):
+        verts = np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0]], dtype=np.float64)
+        faces = np.array([], dtype=np.int32).reshape(0, 3)
+        assert subject.bisect_mesh_plane_vf(verts, faces, plane_z=0.0) == []
+
+    def test_bisect_precision(self):
+        verts, faces = _cube_verts_faces(size=2.0)
+        segments = subject.bisect_mesh_plane_vf(verts, faces, plane_z=0.0, precision=6)
+        for start, end in segments:
+            for coord in start + end:
+                assert round(coord, 6) == coord
+
+
+class TestDissolveFaces:
+    def test_dissolve_cube_into_ngons(self):
+        """A triangulated cube (12 triangles) should dissolve into 6 quad faces."""
+        verts, faces = _cube_verts_faces(size=2.0)
+        edges = np.array(
+            [
+                [0, 1],
+                [1, 2],
+                [2, 3],
+                [3, 0],
+                [4, 5],
+                [5, 6],
+                [6, 7],
+                [7, 4],
+                [0, 4],
+                [1, 5],
+                [2, 6],
+                [3, 7],
+            ],
+            dtype=np.int32,
+        )
+        ngons = subject.dissolve_faces(verts, faces, edges)
+        assert len(ngons) == 6
+        for ngon in ngons:
+            assert len(ngon) == 4
+
+    def test_dissolve_no_edges_returns_triangles(self):
+        """With no original edges, triangles should be returned as-is."""
+        verts, faces = _cube_verts_faces(size=2.0)
+        edges = np.array([], dtype=np.int32).reshape(0, 2)
+        ngons = subject.dissolve_faces(verts, faces, edges)
+        assert len(ngons) == 12
+        for ngon in ngons:
+            assert len(ngon) == 3
+
+    def test_dissolve_empty_faces(self):
+        verts = np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0]], dtype=np.float64)
+        faces = np.array([], dtype=np.int32).reshape(0, 3)
+        edges = np.array([], dtype=np.int32).reshape(0, 2)
+        assert subject.dissolve_faces(verts, faces, edges) == []

src/ifcopenshell-python/test/util/test_space.py

@@ -0,0 +1,186 @@
+# IfcOpenShell - IFC toolkit and geometry engine
+# Copyright (C) 2026 Dion Moult <dion@thinkmoult.com>
+#
+# This file is part of IfcOpenShell.
+#
+# IfcOpenShell is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# IfcOpenShell is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with IfcOpenShell.  If not, see <http://www.gnu.org/licenses/>.
+
+import ifcopenshell.api.geometry
+import ifcopenshell.api.root
+import ifcopenshell.geom
+import ifcopenshell.guid
+import ifcopenshell.util.shape
+import ifcopenshell.util.space as subject
+import pytest
+import shapely
+import test.bootstrap
+
+
+def _build_shapes_dict(ifc_file, elements):
+    """Build a shapes dict as expected by ifcopenshell.util.space functions."""
+    settings = ifcopenshell.geom.settings()
+    settings.set("disable-opening-subtractions", True)
+    settings.set("use-world-coords", True)
+    shapes = {}
+    for element in elements:
+        shape = ifcopenshell.geom.create_shape(settings, element)
+        verts = ifcopenshell.util.shape.get_shape_vertices(shape, shape.geometry)
+        faces = ifcopenshell.util.shape.get_faces(shape.geometry)
+        zs = verts[:, 2]
+        shapes[element.id()] = {
+            "verts": verts,
+            "faces": faces,
+            "bottom_z": float(zs.min()),
+            "top_z": float(zs.max()),
+        }
+    return shapes
+
+
+def _add_extruded_body(ifc_file, element, coords_2d, depth, z_offset=0.0):
+    """Add a body representation (extruded polyline) to an element."""
+    if not ifc_file.by_type("IfcProject"):
+        ifcopenshell.api.root.create_entity(ifc_file, ifc_class="IfcProject")
+    ctx = ifc_file.createIfcGeometricRepresentationContext(
+        ContextType="Model",
+        CoordinateSpaceDimension=3,
+        Precision=1e-5,
+        WorldCoordinateSystem=ifc_file.createIfcAxis2Placement3D(
+            ifc_file.createIfcCartesianPoint((0.0, 0.0, 0.0)),
+            ifc_file.createIfcDirection((0.0, 0.0, 1.0)),
+            ifc_file.createIfcDirection((1.0, 0.0, 0.0)),
+        ),
+    )
+    sub_ctx = ifc_file.createIfcGeometricRepresentationSubContext(
+        ContextIdentifier="Body",
+        ContextType="Model",
+        ParentContext=ctx,
+        TargetView="MODEL_VIEW",
+    )
+    pts = [ifc_file.createIfcCartesianPoint((float(x), float(y))) for x, y in coords_2d]
+    polyline = ifc_file.createIfcPolyline(pts)
+    profile = ifc_file.create_entity("IfcArbitraryClosedProfileDef", ProfileType="CURVE", OuterCurve=polyline)
+    placement = ifc_file.createIfcAxis2Placement3D(
+        ifc_file.createIfcCartesianPoint((0.0, 0.0, z_offset)),
+        ifc_file.createIfcDirection((0.0, 0.0, 1.0)),
+        ifc_file.createIfcDirection((1.0, 0.0, 0.0)),
+    )
+    direction = ifc_file.createIfcDirection((0.0, 0.0, 1.0))
+    solid = ifc_file.createIfcExtrudedAreaSolid(profile, placement, direction, depth)
+    rep = ifc_file.create_entity(
+        "IfcShapeRepresentation",
+        ContextOfItems=sub_ctx,
+        RepresentationIdentifier="Body",
+        RepresentationType="SweptSolid",
+        Items=[solid],
+    )
+    ifcopenshell.api.geometry.assign_representation(ifc_file, product=element, representation=rep)
+
+
+class TestGetBoundaryLines(test.bootstrap.IFC4):
+    def test_returns_segments_for_intersecting_walls(self):
+        wall = ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcWall")
+        _add_extruded_body(self.file, wall, [[-5, -5], [5, -5], [5, 5], [-5, 5]], 3.0)
+        shapes = _build_shapes_dict(self.file, [wall])
+        lines, bounding = subject.get_boundary_lines(self.file, shapes, cut_z=1.0)
+        assert len(lines) > 0
+        assert wall in bounding
+
+    def test_skips_elements_not_intersecting_plane(self):
+        wall = ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcWall")
+        _add_extruded_body(self.file, wall, [[-1, -1], [1, -1], [1, 1], [-1, 1]], 1.0)
+        shapes = _build_shapes_dict(self.file, [wall])
+        lines, bounding = subject.get_boundary_lines(self.file, shapes, cut_z=10.0)
+        assert lines == []
+        assert bounding == []
+
+    def test_skips_non_bounding_classes(self):
+        slab = ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcSlab")
+        _add_extruded_body(self.file, slab, [[-1, -1], [1, -1], [1, 1], [-1, 1]], 1.0)
+        shapes = _build_shapes_dict(self.file, [slab])
+        lines, bounding = subject.get_boundary_lines(self.file, shapes, cut_z=0.5)
+        assert slab not in bounding
+
+
+class TestGetSpacePolygon(test.bootstrap.IFC4):
+    def test_finds_containing_polygon(self):
+        lines = [
+            shapely.LineString([(0, 0), (10, 0)]),
+            shapely.LineString([(10, 0), (10, 10)]),
+            shapely.LineString([(10, 10), (0, 10)]),
+            shapely.LineString([(0, 10), (0, 0)]),
+        ]
+        polygon, _ = subject.get_space_polygon(lines, 5, 5)
+        assert not isinstance(polygon, str)
+        assert polygon.area == pytest.approx(100)
+
+    def test_no_polygons_found(self):
+        polygon, _ = subject.get_space_polygon([], 0, 0)
+        assert polygon == "NO POLYGONS FOUND"
+
+    def test_no_polygon_for_point(self):
+        lines = [
+            shapely.LineString([(0, 0), (10, 0)]),
+            shapely.LineString([(10, 0), (10, 10)]),
+            shapely.LineString([(10, 10), (0, 10)]),
+            shapely.LineString([(0, 10), (0, 0)]),
+        ]
+        polygon, _ = subject.get_space_polygon(lines, 50, 50)
+        assert polygon == "NO POLYGON FOR POINT"
+
+
+class TestGetAutoSpaceHeight(test.bootstrap.IFC4):
+    def test_height_from_top_connection(self):
+        wall = ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcWall")
+        slab = ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcSlab")
+        _add_extruded_body(self.file, wall, [[-5, -5], [5, -5], [5, 5], [-5, 5]], 3.0)
+        _add_extruded_body(self.file, slab, [[-5, -5], [5, -5], [5, 5], [-5, 5]], 0.3, z_offset=3.0)
+        self.file.createIfcRelConnectsElements(
+            GlobalId=ifcopenshell.guid.new(),
+            RelatingElement=slab,
+            RelatedElement=wall,
+            Description="TOP",
+        )
+        shapes = _build_shapes_dict(self.file, [wall, slab])
+        space_polygon = shapely.box(-5, -5, 5, 5)
+        height = subject.get_auto_space_height(self.file, shapes, space_polygon, 0.0, [wall])
+        assert height is not None
+        assert height == pytest.approx(3.0, abs=0.1)
+
+    def test_height_from_elements_above_without_top_connection(self):
+        wall = ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcWall")
+        slab = ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcSlab")
+        _add_extruded_body(self.file, wall, [[-5, -5], [5, -5], [5, 5], [-5, 5]], 3.0)
+        _add_extruded_body(self.file, slab, [[-5, -5], [5, -5], [5, 5], [-5, 5]], 0.3, z_offset=3.0)
+        shapes = _build_shapes_dict(self.file, [wall, slab])
+        space_polygon = shapely.box(-5, -5, 5, 5)
+        height = subject.get_auto_space_height(self.file, shapes, space_polygon, 0.0, [wall])
+        assert height is not None
+        assert height == pytest.approx(3.0, abs=0.1)
+
+    def test_height_from_wall_tops_when_no_slab(self):
+        wall = ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcWall")
+        _add_extruded_body(self.file, wall, [[-5, -5], [5, -5], [5, 5], [-5, 5]], 3.0)
+        shapes = _build_shapes_dict(self.file, [wall])
+        space_polygon = shapely.box(-5, -5, 5, 5)
+        height = subject.get_auto_space_height(self.file, shapes, space_polygon, 0.0, [wall])
+        assert height is not None
+        assert height == pytest.approx(3.0, abs=0.1)
+
+    def test_returns_none_when_no_elements_above(self):
+        wall = ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcWall")
+        _add_extruded_body(self.file, wall, [[-5, -5], [5, -5], [5, 5], [-5, 5]], 3.0)
+        shapes = _build_shapes_dict(self.file, [wall])
+        space_polygon = shapely.box(-100, -100, -90, -90)
+        height = subject.get_auto_space_height(self.file, shapes, space_polygon, 0.0, [])
+        assert height is None