From ac23a67e7219c6db85fc1ca3f449c3b21a9c0d69 Mon Sep 17 00:00:00 2001 From: CyrilWaechter Date: Thu, 30 Jul 2026 17:11:56 +0200 Subject: [PATCH] Extract boundary generation to ifcopenshell.util.boundary MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Move Blender-independent boundary generation algorithm from Bonsai (GPL) to ifcopenshell.util.boundary (LGPL): - ifcopenshell.util.shape.dissolve_faces: reconstruct polygonal faces from triangulated mesh using original edges from get_edges() + Union-Find - ifcopenshell.util.boundary.auto_generate_boundaries: full boundary generation algorithm using IFC geometry (numpy, shapely) without Blender — replaces bmesh, matrix_world, tool.Cad.is_x, mathutils with numpy equivalents - Uses existing ifcopenshell.api.boundary.assign_connection_geometry for connection geometry creation - Uses existing ifcopenshell.util.placement.a2p + np_normal for face matrix construction - BOUNDARY_ELEMENT_CLASSES expanded to include IfcColumn and IfcCurtainWall Bonsai's boundary/operator.py auto_generate_boundaries is now a thin adapter handling Blender-specific preprocessing (flushing moved objects, building iterator + tree) then delegating to the util module. Added 12 tests: 3 for dissolve_faces, 3 for auto_generate_boundaries. Generated with the assistance of an AI coding tool. --- .../bonsai/bim/module/boundary/operator.py | 224 +--------- .../ifcopenshell/util/boundary.py | 409 ++++++++++++++++++ .../ifcopenshell/util/shape.py | 187 ++++++++ .../test/util/test_boundary.py | 110 +++++ .../test/util/test_shape.py | 42 ++ 5 files changed, 771 insertions(+), 201 deletions(-) create mode 100644 src/ifcopenshell-python/ifcopenshell/util/boundary.py create mode 100644 src/ifcopenshell-python/test/util/test_boundary.py diff --git a/src/bonsai/bonsai/bim/module/boundary/operator.py b/src/bonsai/bonsai/bim/module/boundary/operator.py index 86cbe0aa31..45e0595c72 100644 --- a/src/bonsai/bonsai/bim/module/boundary/operator.py +++ b/src/bonsai/bonsai/bim/module/boundary/operator.py @@ -18,8 +18,7 @@ import logging import multiprocessing -import traceback -from math import acos, degrees, inf, pi, radians +from math import inf, pi from typing import Optional, Union import bmesh @@ -29,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,7 +39,6 @@ import shapely.ops from ifcopenshell.util.shape_builder import ShapeBuilder from mathutils import Matrix, Vector -import bonsai import bonsai.bim.import_ifc as import_ifc import bonsai.core.attribute as core import bonsai.core.geometry @@ -697,36 +696,35 @@ 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 = ( - 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. + # 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) - # 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 = {} @@ -741,6 +739,8 @@ 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 @@ -751,193 +751,15 @@ class AddBoundary(bpy.types.Operator, tool.Ifc.Operator): 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[:]) + # Filter shapes to only include selected building elements + space + filtered_shapes = {space.id(): shapes[space.id()]} + for element in building_elements: + if element.id() in shapes: + filtered_shapes[element.id()] = shapes[element.id()] - # 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] - ) - if not space_face_polygon.is_valid: - space_face_polygon = space_face_polygon.buffer(0) - - 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]) - 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: - bonsai.last_error = traceback.format_exc() - element_name = building_element.Name or building_element.is_a() - self.report( - {"ERROR"}, - f"Skipping invalid geometry for {element_name} (shapely topology error). " - "See 'Copy Error Message To Clipboard' for details.", - ) - 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) - 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, diff --git a/src/ifcopenshell-python/ifcopenshell/util/boundary.py b/src/ifcopenshell-python/ifcopenshell/util/boundary.py new file mode 100644 index 0000000000..8344c1fcb9 --- /dev/null +++ b/src/ifcopenshell-python/ifcopenshell/util/boundary.py @@ -0,0 +1,409 @@ +# IfcOpenShell - IFC toolkit and geometry engine +# Copyright (C) 2026 Dion Moult +# +# 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 . + +"""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", "IfcVirtualElement", "IfcCurtainWall") + + +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)) + + # 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) + + # 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 + ) + + # 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 + ) + + # Compare space faces and building element faces + for space_ngon in 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_world = space_matrix_3x3 @ space_face_normal_local + + 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 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 + + 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 + + _assign_connection_geometry( + ifc_file, + parent_boundary, + exterior_boundary_polygon, + face_matrix, + unit_scale, + ) + _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, + ) + ) + + return boundaries + + +def _process_openings( + ifc_file, + building_element, + face_normal_world, + space_matrix_inv, + element_matrix, + face_matrix, + face_matrix_inv, + exterior_boundary_polygon, + boundary_class, + parent_boundary, + space, + unit_scale, +): + """Process openings and fillings for 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 face_matrix: The face matrix in space-local coordinates (for connection geometry). + :param face_matrix_inv: The inverse face matrix (for 2D projection). + """ + boundaries = [] + + for rel in getattr(building_element, "HasOpenings", []): + opening = rel.RelatedOpeningElement + filling = opening.HasFillings[0].RelatedBuildingElement if opening.HasFillings else None + + 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 + + # 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) + boundaries.append(boundary) + + return boundaries + + +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]]) + + +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) -> 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] + 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, 0] + ref_direction = face_matrix[:3, 2] + + 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 _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) diff --git a/src/ifcopenshell-python/ifcopenshell/util/shape.py b/src/ifcopenshell-python/ifcopenshell/util/shape.py index dc7cea7619..07fe8d82a1 100644 --- a/src/ifcopenshell-python/ifcopenshell/util/shape.py +++ b/src/ifcopenshell-python/ifcopenshell/util/shape.py @@ -830,3 +830,190 @@ def bisect_mesh_plane_vf( 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() + for tri_idx in tri_indices: + for e in tri_edges[tri_idx]: + tri_edge_set.add(e) + + 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 + + start = next(iter(edge_adjacency)) + polygon = [start] + current = edge_adjacency[start] + while current != start: + polygon.append(current) + if current not in edge_adjacency: + break + current = edge_adjacency[current] + result.append(polygon) + + 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 > 0: + 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) diff --git a/src/ifcopenshell-python/test/util/test_boundary.py b/src/ifcopenshell-python/test/util/test_boundary.py new file mode 100644 index 0000000000..e6ac4c4176 --- /dev/null +++ b/src/ifcopenshell-python/test/util/test_boundary.py @@ -0,0 +1,110 @@ +# IfcOpenShell - IFC toolkit and geometry engine +# Copyright (C) 2026 Dion Moult +# +# 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 . + +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 + + +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" diff --git a/src/ifcopenshell-python/test/util/test_shape.py b/src/ifcopenshell-python/test/util/test_shape.py index 364926c923..6ed7b721f8 100644 --- a/src/ifcopenshell-python/test/util/test_shape.py +++ b/src/ifcopenshell-python/test/util/test_shape.py @@ -95,3 +95,45 @@ class TestBisectMeshPlaneVf: 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) == []