From cfb7d026d1f4818de7cd5a4c30da8278b77aa751 Mon Sep 17 00:00:00 2001 From: Dion Moult Date: Wed, 26 Feb 2025 18:13:20 +1100 Subject: [PATCH] IfcPatch recipe for Revit 2025 TINs import workaround To understand what's going on, please read the docstring. Then please take as much alcohol required to erase that memory. --- .../ifcpatch/recipes/FixRevit2025TINs.py | 540 ++++++++++++++++++ 1 file changed, 540 insertions(+) create mode 100644 src/ifcpatch/ifcpatch/recipes/FixRevit2025TINs.py diff --git a/src/ifcpatch/ifcpatch/recipes/FixRevit2025TINs.py b/src/ifcpatch/ifcpatch/recipes/FixRevit2025TINs.py new file mode 100644 index 0000000000..f6f375ea23 --- /dev/null +++ b/src/ifcpatch/ifcpatch/recipes/FixRevit2025TINs.py @@ -0,0 +1,540 @@ +# IfcPatch - IFC patching utiliy +# Copyright (C) 2023 Dion Moult +# +# This file is part of IfcPatch. +# +# IfcPatch 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. +# +# IfcPatch 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 IfcPatch. If not, see . + + +import ifcopenshell +import logging +from typing import Optional + + +class Patcher: + def __init__( + self, file: None, logger: logging.Logger, filepath: str, is_solid: bool = True, should_create_edges: bool = True + ): + """Fix missing or spot-coordinate bugged TINs loading in Revit + + TINs exported from 12D or Civil 3D may contain dense or highly obtuse + triangles. Although these will load in Revit, you will not be able to + use Revit's Spot Coordinate or Spot Elevation tool. + + See bug: https://github.com/Autodesk/revit-ifc/issues/511 + + If `is_solid` is enabled, we assume the surface is represented as a + solid (e.g. I have come across surfaces which are extruded by 1mm from + Civil 3D). The solution will delete any faces with a Z normal less + than 0.5. In case the mesh has any side faces or thickness, this should + leave only the top surface which is relevant for spot coordinates and + elevations. + + Vertices closer than 10mm will also be merged to prevent dense + portions of the TIN at a minor sacrifice of surveying accuracy. It will + also triangulate all meshes to prevent non-coplanar surfaces, and + delete any obtuse triangles where one of their XY angles is less than + 0.3 degrees. Therefore the result will contain some minor "holes" in + the TIN, but these holes will only be in dense triangles that Revit + can't handle anyway and won't affect most coordination tasks. + + After that, it will: + + 1. Reassign everything to an IfcGeographicElement + 2. Detect boundary edges and create an edge-only IfcVirtualElement. + Good for clean viz in Revit. + 3. Create a copy of the object which has no sharp faces. This will + allow Revit's spot coordinate tool to work on any arbitrary face + surface. Note that Revit cannot snap to edges or vertices on this + object. + 4. Create a copy of the obejct which has one artifically injected sharp + face. This trick allows Revit's spot coordinate tool to snap to + edges and points. However, Revit cannot sample an arbitrary + surface. By combining this object with the previous object, you get + the best of both worlds. + + If you're thinking that this overlapping, Z-fighting, duplication of + objects with arbitrary almost-degenerate triangles being added is a + horrific abomination in the world of software workarounds, you are + absolutely correct. + + This is a variation of FixRevitTINs which has been tested on Revit <= + 2023. I've tested this one on Revit 2025 (the behaviour has changed). + + Note that you may may want to run other tools like + OffsetObjectPlacements or ResetAbsoluteCoordinates to fix large + coordinates as these can also cause issues in Revit (such as inaccuracy + or inability to use the Spot Coordinate / Elevation tool). + + This patch is designed to work on any TIN-like export, typically coming + from civil software. It also requires you to run it using Blender, as + the geometric modification uses the Blender geometry engine. + + `filepath` argument is required for this recipe, `file` argument is + ignored. + + :param filepath: The filepath of the IFC model. This is required to + load into Bonsai. + :param is_solid: If true, assume a thickness and delete anything that + isn't the top face. + :param should_create_edges: If true, a new IfcVirtualElement is created + representing the perimeter of the objects. This allows you to to + hide regular surface edges in Revit and only use the perimeter edge + for visualisation. + :filter_glob filepath: *.ifc;*.ifczip;*.ifcxml + + Example: + + .. code:: python + + ifcpatch.execute({"input": "input.ifc", "recipe": "FixRevit2025TINs", "arguments": []}) + """ + self.file = file + self.filepath = filepath + self.logger = logger + self.is_solid = is_solid + self.should_create_edges = should_create_edges + + def patch(self) -> None: + import bpy + import bmesh + import bonsai.tool as tool + import ifcopenshell.util.shape_builder + + bpy.context.scene.BIMProjectProperties.should_use_native_meshes = True + bpy.ops.bim.load_project(filepath=self.filepath) + + old_history_size = tool.Ifc.get().history_size + old_undo_steps = bpy.context.preferences.edit.undo_steps + tool.Ifc.get().history_size = 0 + bpy.context.preferences.edit.undo_steps = 0 + + self.unit_scale = ifcopenshell.util.unit.calculate_unit_scale(tool.Ifc.get()) + + for obj in bpy.data.objects: + if not obj.BIMObjectProperties.ifc_definition_id or not obj.data: + continue + if not obj.data.polygons: + continue + element = tool.Ifc.get_entity(obj) + element.PredefinedType = "USERDEFINED" + element.ObjectType = "TIN" + element = ifcopenshell.util.schema.reassign_class(tool.Ifc.get(), element, "IfcGeographicElement") + + bm = bmesh.new() + bm.from_mesh(obj.data) + faces_to_delete = [] + if self.is_solid: + for face in bm.faces: + global_normal = obj.matrix_world.to_3x3() @ face.normal + if global_normal.z < 0.5: + faces_to_delete.append(face) + bmesh.ops.delete(bm, geom=faces_to_delete, context="FACES_ONLY") + bmesh.ops.triangulate(bm, faces=bm.faces[:], quad_method="BEAUTY", ngon_method="BEAUTY") + bm.to_mesh(obj.data) + bm.free() + obj.data.update() + + if self.should_create_edges: + self.create_edges(obj) + + self.create_face_sampleable_object(obj) + self.create_edge_sampleable_object(obj) + + tool.Ifc.get().history_size = old_history_size + bpy.context.preferences.edit.undo_steps = old_undo_steps + + self.file = tool.Ifc.get() + + def create_edges(self, obj): + import bpy + import bmesh + import bonsai.tool as tool + import ifcopenshell.util.element + import ifcopenshell.util.representation + import ifcopenshell.api.root + import ifcopenshell.api.type + import ifcopenshell.api.spatial + import ifcopenshell.api.geometry + + element = tool.Ifc.get_entity(obj) + data = obj.data + bm = bmesh.new() + bm.from_mesh(data) + + bm.faces.ensure_lookup_table() + + if self.is_solid: + faces_to_delete = [] + for face in bm.faces: + if face.normal.z < 0.5: + faces_to_delete.append(face) + bmesh.ops.delete(bm, geom=faces_to_delete, context="FACES_ONLY") + + bmesh.ops.remove_doubles(bm, verts=bm.verts, dist=0.01) + + edges_to_delete = [] + bm.faces.ensure_lookup_table() + for edge in bm.edges: + if len(edge.link_faces) != 1: + edges_to_delete.append(edge) + + bmesh.ops.delete(bm, geom=edges_to_delete, context="EDGES_FACES") + bmesh.ops.remove_doubles(bm, verts=bm.verts, dist=0.0001) + + bm.verts.ensure_lookup_table() + for vert in bm.verts: + vert.co.z += 0.003 + + mesh = bpy.data.meshes.new("Mesh") + bm.to_mesh(mesh) + + obj = bpy.data.objects.new("Perimeter", mesh) + bpy.context.scene.collection.objects.link(obj) + with bpy.context.temp_override(**tool.Blender.get_viewport_context()): + tool.Blender.select_and_activate_single_object(bpy.context, obj) + bpy.ops.object.convert(target="CURVE") + + context = ifcopenshell.util.representation.get_context(tool.Ifc.get(), "Model", "Body", "MODEL_VIEW") + + if self.file.schema == "IFC2X3": + curves = self.create_curves_from_curve_ifc2x3(is_2d=False, curve_object_data=obj.data) + else: + curves = self.create_curves_from_curve(is_2d=False, curve_object_data=obj.data) + + representation = tool.Ifc.get().createIfcShapeRepresentation( + context, context.ContextIdentifier, "Curve3D", curves + ) + + element2 = ifcopenshell.api.root.copy_class(tool.Ifc.get(), product=element) + element2.Name += "-boundary" + element2.ObjectType = "TINBOUNDARY" + ifcopenshell.api.geometry.assign_representation(tool.Ifc.get(), element2, representation) + ifcopenshell.util.schema.reassign_class(tool.Ifc.get(), element2, "IfcVirtualElement") + bm.free() + + def create_face_sampleable_object(self, obj): + # No sharp faces + import bpy + import bmesh + import bonsai.tool as tool + import ifcopenshell.util.element + import ifcopenshell.util.representation + import ifcopenshell.api.root + import ifcopenshell.api.type + import ifcopenshell.api.spatial + import ifcopenshell.api.geometry + from math import degrees + + print("working on ", obj.name) + element = tool.Ifc.get_entity(obj) + data = obj.data + bm = bmesh.new() + bm.from_mesh(data) + + bm.faces.ensure_lookup_table() + + angle_threshold = 0.3 + for polygon in bm.faces: + try: + v1, v2, v3 = [v.co for v in polygon.verts] + d1 = degrees((v2 - v1).angle(v3 - v1)) + d2 = degrees((v3 - v2).angle(v1 - v2)) + d3 = degrees((v1 - v3).angle(v2 - v3)) + if d1 < angle_threshold or d2 < angle_threshold or d3 < angle_threshold: + print("removing", d1, d2, d3, polygon) + bm.faces.remove(polygon) + except: + print("removing", polygon) + bm.faces.remove(polygon) + + context = ifcopenshell.util.representation.get_context(tool.Ifc.get(), "Model", "Body", "MODEL_VIEW") + + builder = ifcopenshell.util.shape_builder.ShapeBuilder(tool.Ifc.get()) + verts = [v.co / self.unit_scale for v in bm.verts] + faces = [[v.index for v in p.verts] for p in bm.faces] + item = builder.mesh(verts, faces) + representation = builder.get_representation(context, [item]) + element2 = ifcopenshell.api.root.copy_class(tool.Ifc.get(), product=element) + element2.Name += "-face-sample" + print("new", element2, element) + ifcopenshell.api.geometry.assign_representation(tool.Ifc.get(), element2, representation) + + bm.free() + + def create_edge_sampleable_object(self, obj): + # This is crazy but we need a sharp face per island + import bpy + import bmesh + import bonsai.tool as tool + import ifcopenshell.util.element + import ifcopenshell.util.representation + import ifcopenshell.api.root + import ifcopenshell.api.type + import ifcopenshell.api.spatial + import ifcopenshell.api.geometry + from math import degrees, radians, sin + from mathutils import Matrix + + # Get the active object (assumed to have a mesh) + mesh = obj.data + + # Create a BMesh representation + bm = bmesh.new() + bm.from_mesh(mesh) + + # First purge existing sharp edges (don't ask me why, really don't) + bm.faces.ensure_lookup_table() + angle_threshold = 0.3 + for polygon in bm.faces: + try: + # v1, v2, v3 = [v.co.to_2d() for v in polygon.verts] + v1, v2, v3 = [v.co for v in polygon.verts] + d1 = degrees((v2 - v1).angle(v3 - v1)) + d2 = degrees((v3 - v2).angle(v1 - v2)) + d3 = degrees((v1 - v3).angle(v2 - v3)) + if d1 < angle_threshold or d2 < angle_threshold or d3 < angle_threshold: + print("removing", d1, d2, d3, polygon) + bm.faces.remove(polygon) + except: + print("removing", polygon) + bm.faces.remove(polygon) + bm.faces.ensure_lookup_table() + + # Now we add our own. + + # A set to mark faces that have been visited + visited_faces = set() + + def get_island(start_face): + """Return the connected set of faces (a 'mesh island') starting from start_face.""" + island = set() + stack = [start_face] + while stack: + f = stack.pop() + if f in island: + continue + island.add(f) + for edge in f.edges: + # For every face sharing this edge, add to the stack + for f2 in edge.link_faces: + if f2 not in island: + stack.append(f2) + return island + + # Loop over all faces and process each island once + islands_count = 0 + for face in bm.faces: + if face in visited_faces: + continue + + # Get the connected component (island) containing this face + island = get_island(face) + visited_faces |= island # mark all island faces as visited + islands_count += 1 + + boundary_edge = None + for face in island: + for edge in face.edges: + # Count how many faces in the island use this edge. + count = sum(1 for f in edge.link_faces if f in island) + if count == 1: + boundary_edge = edge + break + if boundary_edge: + break + + # If no boundary edge is found, skip this island. + if boundary_edge is None: + continue + + # Use the two vertices of the boundary edge as A and B. + A, B = boundary_edge.verts[0], boundary_edge.verts[1] + + # THIRD ATTEMPT + + # Use the normal from the boundary face (i.e. the single linked face of the boundary edge) + base_face = boundary_edge.link_faces[0] + plane_normal = ( + base_face.normal.copy() + ) # This normal defines the plane in which we'll construct the triangle + + # Compute the edge AB vector and its length. + AB_vec = B.co - A.co + d = AB_vec.length + if d == 0: + continue # degenerate edge + + # --- Desired angles for the new triangle (in degrees) --- + # This is the crazy degenerate triangle + angle_A_deg = 0.004 # angle at vertex A + angle_B_deg = 1.146 # angle at vertex B + angle_C_deg = 178.85 # angle at vertex C (note: 180 - (0.004 + 1.146) = 178.85) + + # Convert angles to radians. + angle_A = radians(angle_A_deg) + angle_B = radians(angle_B_deg) + angle_C = radians(angle_C_deg) + + # --- Use the law of sines to compute the new triangle's side lengths --- + # In triangle ABC, with AB opposite angle C: + # AB / sin(angle_C) = AC / sin(angle_B) = BC / sin(angle_A) + # We'll compute AC (from A) as it is needed to place the new vertex. + AC_length = d * sin(angle_B) / sin(angle_C) + + # --- Determine the direction for AC in the triangle's plane --- + # Starting at A, the direction of AB is our baseline. + u = AB_vec.normalized() + + # To get the direction for AC, rotate u by angle_A about the plane_normal. + rot_mat = Matrix.Rotation(angle_A, 3, plane_normal) + dA = u.copy() + dA.rotate(rot_mat) + + # Compute the position for the new vertex C. + C_co = A.co + AC_length * dA + + # --- Create the new vertex and triangle face in the BMesh --- + new_vert = bm.verts.new(C_co) + bm.verts.index_update() # update indices if needed + + # Create the new triangle face from vertices A, B, and new_vert. + # (The order of vertices may be adjusted if you need a specific winding.) + try: + new_face = bm.faces.new((A, B, new_vert)) + visited_faces.add(new_face) + print("added new face") + except ValueError: + # Face already exists or some error occurred + print("Could not create face on island", islands_count) + + context = ifcopenshell.util.representation.get_context(tool.Ifc.get(), "Model", "Body", "MODEL_VIEW") + builder = ifcopenshell.util.shape_builder.ShapeBuilder(tool.Ifc.get()) + verts = [v.co / self.unit_scale for v in bm.verts] + faces = [[v.index for v in p.verts] for p in bm.faces] + item = builder.mesh(verts, faces) + representation = builder.get_representation(context, [item]) + element = tool.Ifc.get_entity(obj) + element2 = ifcopenshell.api.root.copy_class(tool.Ifc.get(), product=element) + element2.Name += "-edge-sample" + print("new", element2, element) + ifcopenshell.api.geometry.assign_representation(tool.Ifc.get(), element2, representation) + + bm.free() + + print("Added a triangle to", islands_count, "mesh island(s).") + return + + def create_curves_from_curve_ifc2x3( + self, is_2d: bool = False, curve_object_data=None + ) -> list[ifcopenshell.entity_instance]: + import bonsai.tool as tool + + dim = (lambda v: v.xy) if is_2d else (lambda v: v.xyz) + results = [] + for spline in curve_object_data.splines: + points = spline.bezier_points[:] + spline.points[:] + if spline.use_cyclic_u: + points.append(points[0]) + ifc_points = [self.create_cartesian_point(*dim(point.co)) for point in points] + results.append(tool.Ifc.get().createIfcPolyline(ifc_points)) + return results + + def create_curves_from_curve( + self, is_2d: bool = False, curve_object_data=None + ) -> list[ifcopenshell.entity_instance]: + import bonsai.tool as tool + import numpy as np + + dim = (lambda v: v.xy) if is_2d else (lambda v: v.xyz) + to_units = lambda v: np.array([self.convert_si_to_unit(i) for i in v]) + builder = ifcopenshell.util.shape_builder.ShapeBuilder(tool.Ifc.get()) + results = [] + + for spline in curve_object_data.splines: + points = spline.bezier_points[:] + spline.points[:] + + points = [to_units(dim(p.co)) for p in points] + closed_polyline = spline.use_cyclic_u and len(points) > 1 + results.append(builder.polyline(points, closed=closed_polyline)) + + return results + + def create_cartesian_point( + self, x: float, y: float, z: Optional[float] = None, is_model_coords: bool = True + ) -> ifcopenshell.entity_instance: + """Create IfcCartesianPoint. + + x, y, z coords are provided in SI units. + """ + x = self.convert_si_to_unit(x) + y = self.convert_si_to_unit(y) + z = self.convert_si_to_unit(z) + return self.file.createIfcCartesianPoint((x, y, z)) + + def create_curves_from_mesh(self, geom_data) -> list[ifcopenshell.entity_instance]: + curves = [] + points = self.create_cartesian_point_list_from_vertices(geom_data.vertices) + edge_loops = [] + previous_edge = None + edge_loop = [] + for i, edge in enumerate(geom_data.edges): + if previous_edge is None: + edge_loop = [self.file.createIfcLineIndex((edge.vertices[0] + 1, edge.vertices[1] + 1))] + elif edge.vertices[0] == previous_edge.vertices[1]: + edge_loop.append(self.file.createIfcLineIndex((edge.vertices[0] + 1, edge.vertices[1] + 1))) + else: + edge_loops.append(edge_loop) + edge_loop = [self.file.createIfcLineIndex((edge.vertices[0] + 1, edge.vertices[1] + 1))] + previous_edge = edge + edge_loops.append(edge_loop) + for edge_loop in edge_loops: + curves.append(self.file.createIfcIndexedPolyCurve(points, edge_loop)) + return curves + + def create_curves_from_mesh_ifc2x3(self, geom_data) -> list[ifcopenshell.entity_instance]: + curves = [] + points = [self.create_cartesian_point(v.co.x, v.co.y, v.co.z) for v in geom_data.vertices] + edge_loops = [] + previous_edge = None + edge_loop = [] + for i, edge in enumerate(geom_data.edges): + if previous_edge is None: + edge_loop = [edge.vertices] + elif edge.vertices[0] == previous_edge.vertices[1]: + edge_loop.append(edge.vertices) + else: + edge_loops.append(edge_loop) + edge_loop = [edge.vertices] + previous_edge = edge + edge_loops.append(edge_loop) + for edge_loop in edge_loops: + loop_points = [points[p[0]] for p in edge_loop] + loop_points.append(points[edge_loop[-1][1]]) + curves.append(self.file.createIfcPolyline(loop_points)) + return curves + + def create_cartesian_point_list_from_vertices(self, vertices) -> ifcopenshell.entity_instance: + import numpy as np + from ifcopenshell.util.shape_builder import ifc_safe_vector_type + + # Catch values as floats to benefit from fast buffer copy. + coords = np.empty(len(vertices) * 3, dtype="f") + vertices.foreach_get("co", coords) + coords = coords.reshape(-1, 3) + coords_class = "IfcCartesianPointList3D" + return self.file.create_entity(coords_class, ifc_safe_vector_type(self.convert_si_to_unit(coords))) + + def convert_si_to_unit(self, co): + return co / self.unit_scale