Fix space boundary generation regressions

When several elements match the same space face, offset matches that only
duplicate coplanar coverage are now skipped, and a single bounding element
within a small plane offset gets the full space face instead of a clipped
polygon. Existing boundaries are removed before regeneration so stale 2nd
level boundaries are not left behind, and the Bonsai operator delegates
element filtering to auto_generate_boundaries.

Regenerates SmallHouse boundaries to match the reference output and keeps
the ExternalEarth opening unioning intact.

Generated with the assistance of an AI coding tool.

.gitmodules

@@ -20,3 +20,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/module/boundary/operator.py

@@ -719,11 +719,6 @@ class AddBoundary(bpy.types.Operator, tool.Ifc.Operator):
         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 for elements that already have boundaries
-        for boundary in space.BoundedBy:
-            if boundary.RelatedBuildingElement in building_elements:
-                building_elements.remove(boundary.RelatedBuildingElement)
-
         # Build shapes dict with iterator (parallel, includes space + building elements)
         include = building_elements + [space]
         tree = ifcopenshell.geom.tree()
@@ -745,13 +740,9 @@ class AddBoundary(bpy.types.Operator, tool.Ifc.Operator):
                 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]
-
-        if not building_elements:
-            return "No building elements found to create boundaries."
-
-        # Filter shapes to only include selected building elements + 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:

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

@@ -43,7 +43,27 @@ import shapely.ops
 
 logger = logging.getLogger("ImportIFC")
 
-BOUNDARY_ELEMENT_CLASSES = ("IfcWall", "IfcColumn", "IfcSlab", "IfcVirtualElement", "IfcCurtainWall")
+BOUNDARY_ELEMENT_CLASSES = (
+    "IfcWall",
+    "IfcColumn",
+    "IfcSlab",
+    "IfcRoof",
+    "IfcVirtualElement",
+    "IfcCurtainWall",
+    "IfcWindow",
+    "IfcDoor",
+)
+
+# Distance (in meters) below which an element face is considered coplanar with
+# the space face it bounds. Matches beyond this distance within the larger
+# tolerance are only kept when no coplanar element already covers the face.
+COPLANAR_TOL = 0.05
+
+# 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.2
 
 
 def auto_generate_boundaries(
@@ -78,10 +98,12 @@ def auto_generate_boundaries(
     for ifc_class in boundary_element_classes:
         building_elements.extend(ifc_file.by_type(ifc_class))
 
-    # Don't generate boundaries for elements that already have boundaries
-    for boundary in space.BoundedBy:
+    # 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:
-            building_elements.remove(boundary.RelatedBuildingElement)
+            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]
@@ -107,203 +129,480 @@ def auto_generate_boundaries(
             es["verts"], es["faces"], es["edges"], merge_coplanar=True
         )
 
-    # Compare space faces and building element faces
+    # 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()
-    for space_ngon in space_ngons:
+    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]
-        # Normal from local verts, then transform to world via space placement
         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_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:
-            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)
-
-            for ngon in element_ngons[element.id()]:
-                elem_verts_l = element_shape["verts"][ngon]
-                # Normal from local verts, transform to world via element placement
-                elem_face_normal_local = _face_normal(elem_verts_l)
-                if elem_face_normal_local is None:
-                    continue
-                elem_face_normal_world = element_matrix_3x3 @ elem_face_normal_local
-
-                # Both normals point outward from their respective solids.
-                # Adjacent faces have anti-parallel normals (angle ≈ 180°).
-                # Virtual elements use parallel normals (angle ≈ 0°).
-                angle = degrees(acos(max(min(float(np.dot(space_face_normal_world, elem_face_normal_world)), 1), -1)))
-                if _is_x(angle, 180, tolerance=2):
-                    pass
-                elif element.is_a("IfcVirtualElement") and _is_x(angle, 0, tolerance=2):
-                    pass
-                else:
-                    continue
-
-                # Distance check: transform space vert to element-local, compare to element face
-                # space-local -> world -> element-local
-                space_vert_in_elem = sb.np_apply_matrix(space_verts_l[:1], element_matrix_inv @ space_matrix)[0]
-                dist = float(np.dot(space_vert_in_elem - elem_verts_l[0], elem_face_normal_local))
-                if abs(dist) > 0.05:
-                    continue
-
-                # Build face matrix in space-local coordinates
-                # (assign_connection_geometry expects location/axes relative to space placement)
-                face_matrix = _face_matrix_from_verts(space_verts_l[:3])
-                face_matrix_inv = np.linalg.inv(face_matrix)
-
-                # Project space face (already space-local) to 2D
-                space_face_polygon = _verts_to_polygon(space_verts_l, face_matrix_inv)
-                if not space_face_polygon.is_valid:
-                    space_face_polygon = space_face_polygon.buffer(0)
-
-                # Transform element verts to space-local, then project to 2D
-                # element-local -> world -> space-local
-                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)
-                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 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):
+            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). When the
+        # face is already covered coplanarly, offset matches that only duplicate
+        # that coverage are skipped.
+        coplanar_union = None
+        for _, dist_min, _, gross_boundary_polygon, _ in candidates:
+            if dist_min <= COPLANAR_TOL:
+                coplanar_union = (
+                    gross_boundary_polygon if coplanar_union is None else coplanar_union.union(gross_boundary_polygon)
+                )
+
+        surviving_candidates = []
+        for element, dist_min, plane_offset_min, gross_boundary_polygon, matched_elem_normal in candidates:
+            if dist_min > COPLANAR_TOL and coplanar_union is not None:
+                if gross_boundary_polygon.difference(coplanar_union).area < 1e-4:
                     continue
-                if gross_boundary_polygon.is_empty:
+            surviving_candidates.append(
+                (element, dist_min, plane_offset_min, gross_boundary_polygon, matched_elem_normal)
+            )
+
+        # 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:
+            if dist_min > COPLANAR_TOL and coplanar_union is not None:
+                if gross_boundary_polygon.difference(coplanar_union).area < 1e-4:
                     continue
 
-                exterior_boundary_polygon = shapely.Polygon(gross_boundary_polygon.exterior.coords)
+            exterior_boundary_polygon = shapely.Polygon(gross_boundary_polygon.exterior.coords)
 
-                # Create parent boundary
-                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
+            opening_source_element = element
+            for rel in getattr(element, "Decomposes", []):
+                if rel.RelatingObject.is_a() in BOUNDARY_ELEMENT_CLASSES:
+                    element = rel.RelatingObject
+                    break
 
-                _assign_connection_geometry(
+            # 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,
-                    parent_boundary,
-                    exterior_boundary_polygon,
+                    openings_to_process,
                     face_matrix,
+                    boundary_class,
+                    parent_boundary,
+                    space,
                     unit_scale,
+                    processed_fillings,
+                    covered_by_face[space_ngon_idx],
                 )
-                _set_boundary_name(parent_boundary)
-                boundaries.append(parent_boundary)
+            )
 
-                # Process openings
-                boundaries.extend(
-                    _process_openings(
-                        ifc_file,
-                        element,
-                        elem_face_normal_world,
-                        space_matrix_inv,
-                        element_matrix,
-                        face_matrix,
-                        face_matrix_inv,
-                        exterior_boundary_polygon,
-                        boundary_class,
-                        parent_boundary,
-                        space,
-                        unit_scale,
-                        processed_fillings,
-                    )
-                )
+    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,
-    building_element,
-    face_normal_world,
-    space_matrix_inv,
-    element_matrix,
+    openings_to_process,
     face_matrix,
-    face_matrix_inv,
-    exterior_boundary_polygon,
     boundary_class,
     parent_boundary,
     space,
     unit_scale,
     processed_fillings: set[int],
+    covered_polygons: list,
 ):
-    """Process openings and fillings for a building element.
+    """Create boundaries for the fillings of openings in a building element.
 
-    :param face_normal_world: The building element face normal in world space.
-    :param space_matrix_inv: Inverse of the space placement matrix.
-    :param element_matrix: The building element placement matrix.
+    :param openings_to_process: Tuples of (opening, filling, opening polygon).
     :param face_matrix: The face matrix in space-local coordinates (for connection geometry).
-    :param face_matrix_inv: The inverse face matrix (for 2D projection).
     :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 rel in getattr(building_element, "HasOpenings", []):
-        opening = rel.RelatedOpeningElement
-        filling = opening.HasFillings[0].RelatedBuildingElement if opening.HasFillings else None
-        filling_id = (filling or opening).id()
+    for opening, filling, opening_polygon in openings_to_process:
+        filling_id = filling.id()
         if filling_id in processed_fillings:
             continue
 
-        settings = ifcopenshell.geom.settings()
-        try:
-            shape = ifcopenshell.geom.create_shape(settings, opening)
-        except Exception:
-            continue
-        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]
-            # Normal from local verts, transform to world via opening placement
-            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
-            # Transform opening verts to space-local: opening-local -> world -> space-local
-            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:
-            continue
-
-        opening_polygon = shapely.ops.unary_union(opening_polygons)
-
-        if opening_polygon.intersection(exterior_boundary_polygon).area == 0:
-            continue
-
         boundary = ifcopenshell.api.root.create_entity(ifc_file, ifc_class=boundary_class)
         boundary.RelatingSpace = space
         boundary.RelatedBuildingElement = filling or opening
@@ -325,16 +624,198 @@ def _process_openings(
             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
-    return sb.np_normal([verts[0], verts[1], verts[2]])
+    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:
@@ -345,9 +826,11 @@ def _face_matrix_from_verts(verts3: np.ndarray) -> np.ndarray:
     return ifcopenshell.util.placement.a2p(o=p1, z=z, x=x)
 
 
-def _verts_to_polygon(verts: np.ndarray, face_matrix_inv: np.ndarray) -> shapely.Polygon:
+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])
 
 
@@ -412,6 +895,20 @@ def _set_boundary_name(boundary: ifcopenshell.entity_instance) -> None:
         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/shape.py

@@ -991,7 +991,7 @@ def _merge_coplanar_ngons(
         edge2 = v2 - v0
         normal = np.cross(edge1, edge2)
         norm = np.linalg.norm(normal)
-        if norm > 0:
+        if norm > 1e-8:
             normal = normal / norm
         ngon_normals[root] = normal
 

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

@@ -16,6 +16,13 @@
 # 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
@@ -80,6 +87,126 @@ def _build_shapes_dict(ifc_file, elements):
     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 _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")
@@ -108,3 +235,88 @@ class TestAutoGenerateBoundaries(test.bootstrap.IFC4):
         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]