From 2d3295ec54a800fc4f2055e37ef9303fe5404409 Mon Sep 17 00:00:00 2001 From: Andrej730 Date: Mon, 21 Aug 2023 17:22:20 +0500 Subject: [PATCH] more work on transitions for MEP - support transitions from and to circle profiles - reworked transition length algorithm now it should be more accurate - added support for creating transitions between profiles that are parallel but not collinear --- .../blenderbim/bim/module/model/mep.py | 85 ++++- src/blenderbim/blenderbim/tool/cad.py | 18 +- src/blenderbim/test/tool/test_cad.py | 58 ++++ .../ifcopenshell/util/shape_builder.py | 324 +++++++++++++++--- 4 files changed, 405 insertions(+), 80 deletions(-) create mode 100644 src/blenderbim/test/tool/test_cad.py diff --git a/src/blenderbim/blenderbim/bim/module/model/mep.py b/src/blenderbim/blenderbim/bim/module/model/mep.py index d68c6d43f6..a0282a9a24 100644 --- a/src/blenderbim/blenderbim/bim/module/model/mep.py +++ b/src/blenderbim/blenderbim/bim/module/model/mep.py @@ -29,6 +29,7 @@ import ifcopenshell.util.system import ifcopenshell.util.element import ifcopenshell.util.representation import mathutils.geometry +import numpy as np import blenderbim.bim.handler import blenderbim.core.type import blenderbim.core.root @@ -347,6 +348,8 @@ class MEPGenerator: There lies the problem that it won't be able to identify the fittings that were not yet connected to any segments yet. """ + + # TODO: check angle, start, end and offset for transitions if not isinstance(segment_or_segments, collections.abc.Iterable): segments = [segment_or_segments] ports = [port_or_ports] @@ -436,7 +439,7 @@ class MEPGenerator: if element_type is None: skipped_the_occurrence = True break - + fitting_data.append((element_type, port.PredefinedType, port.SystemType)) # if we skipped the occurrence we still can other occurrences @@ -624,10 +627,13 @@ class MEPAddTransition(bpy.types.Operator, tool.Ifc.Operator): start_axis = tool.Model.get_flow_segment_axis(start_object) end_axis = tool.Model.get_flow_segment_axis(end_object) + start_object_rotation = start_object.matrix_world.to_quaternion() + start_object_z_basis = start_object_rotation.to_matrix().col[2] # z basis vector + keep_only_z_axis = lambda p_ws: p_ws.dot(start_object_z_basis) * start_object_z_basis # TODO: support cases when segments are partially or completely overlapping each other - if not tool.Cad.are_edges_collinear(start_axis, end_axis): - self.report({"ERROR"}, f"Failed to add transition - non collinear segments are not yet supported.") + if not tool.Cad.are_edges_parallel(start_axis, end_axis): + self.report({"ERROR"}, f"Failed to add transition - segments are not parallel.") return {"CANCELLED"} start_segment_data = MEPGenerator().get_segment_data(start_element) @@ -647,42 +653,78 @@ class MEPAddTransition(bpy.types.Operator, tool.Ifc.Operator): (start_segment_data["start_point"], start_segment_data["end_point"]), (end_segment_data["start_point"], end_segment_data["end_point"]), ) - + + # figure profile offset + base_transition_dir = keep_only_z_axis(end_point - start_point).normalized() + flip_profile_offset = base_transition_dir.dot(start_object_z_basis) < 0 + + if tool.Cad.are_edges_collinear(start_axis, end_axis): + profile_offset = None + else: + to_start_object_space = start_object_rotation.inverted() + profile_offset = ( + (to_start_object_space @ end_object.location) - (to_start_object_space @ start_object.location) + ).xy + if tool.Cad.is_x(profile_offset.length_squared, 0): + profile_offset = None + else: + profile_offset = profile_offset / si_conversion + if flip_profile_offset: + profile_offset *= V(1, -1) + + # world space profile offset + profile_offset_ws = ( + start_object_rotation @ (profile_offset * si_conversion).to_3d() if profile_offset else V(0, 0, 0) + ) + + # will need entire_length to check that transition length fill fit first_segment_start, second_segment_end = [ - p for p in ( - start_segment_data["start_point"], - start_segment_data["end_point"], - end_segment_data["start_point"], - end_segment_data["end_point"]) + p + for p in ( + start_segment_data["start_point"], + start_segment_data["end_point"], + end_segment_data["start_point"], + end_segment_data["end_point"], + ) if p not in (start_point, end_point) ] entire_length = (first_segment_start - second_segment_end).length - transition_dir = (end_point - start_point).normalized() + # can't rely on (end_point-start_point) here because + # transition might change the segments length and therefore direction will be changed + segments_dir = (start_point - first_segment_start).normalized() start_port = points_ports_map[start_point] end_port = points_ports_map[end_point] # add transition representation builder = ShapeBuilder(ifc_file) rep, transition_data = builder.mep_transition_shape( - start_element, end_element, self.start_length / si_conversion, self.end_length / si_conversion + start_element, + end_element, + self.start_length / si_conversion, + self.end_length / si_conversion, + profile_offset=profile_offset, ) if not rep: self.report({"ERROR"}, f"Failed to add transition - this kind of profiles is not yet supported.") return {"CANCELLED"} - - # TODO: test it + full_transition_length = transition_data["full_transition_length"] * si_conversion if full_transition_length >= entire_length: - self.report({"ERROR"}, f"Failed to add transition - transition length is larger the segments and the distance between them.") - # TODO: handle the case without creating representation in the first place? + self.report( + {"ERROR"}, + f"Failed to add transition - transition length is larger the segments and the distance between them.\n" + + f"Transition length: {full_transition_length:.2f}m, segments length: {entire_length:.2f}m", + ) + # TODO: handle the case without creating a representation in the first place? ifcopenshell.api.run("geometry.remove_representation", ifc_file, representation=rep) return {"CANCELLED"} - middle_point = (start_point + end_point) / 2 - start_segment_extend_point = middle_point - transition_dir * full_transition_length / 2 - end_segment_extend_point = middle_point + transition_dir * full_transition_length / 2 + middle_point = keep_only_z_axis((start_point + end_point) / 2 - start_point) + start_point + start_segment_extend_point = middle_point - segments_dir * full_transition_length / 2 + end_segment_extend_point = middle_point + segments_dir * full_transition_length / 2 + profile_offset_ws + transition_dir = keep_only_z_axis(end_segment_extend_point - start_segment_extend_point).normalized() DumbProfileJoiner().join_E(start_object, start_segment_extend_point) DumbProfileJoiner().join_E(end_object, end_segment_extend_point) @@ -691,6 +733,10 @@ class MEPAddTransition(bpy.types.Operator, tool.Ifc.Operator): ) transition_type = fitting_data["fitting_type"] if fitting_data else None + if transition_type: + # TODO: handle the case without creating a representation in the first place? + ifcopenshell.api.run("geometry.remove_representation", ifc_file, representation=rep) + start_port_match = fitting_data["start_port_match"] if fitting_data else True if not transition_type: @@ -722,8 +768,9 @@ class MEPAddTransition(bpy.types.Operator, tool.Ifc.Operator): # adjust transition segment rotation and location transition_obj.matrix_world = start_object.matrix_world context.view_layer.update() + transition_obj_dir = tool.Cad.get_edge_direction(tool.Model.get_flow_segment_axis(transition_obj)) - direction_match = tool.Cad.are_vectors_equal(transition_obj_dir, transition_dir) + direction_match = tool.Cad.are_vectors_equal(transition_dir, transition_obj_dir) # if there are no mismatches or everything matches up we don't need to flip the transition if start_port_match != direction_match: diff --git a/src/blenderbim/blenderbim/tool/cad.py b/src/blenderbim/blenderbim/tool/cad.py index 7d0fe1bca7..3443e6cca7 100644 --- a/src/blenderbim/blenderbim/tool/cad.py +++ b/src/blenderbim/blenderbim/tool/cad.py @@ -236,20 +236,16 @@ class Cad: return (edge[1] - edge[0]).normalized() @classmethod - def are_edges_collinear(cls, edge1, edge2): - def is_point_on_line(p, edge): - a1, a2 = edge - # comparing slopes between PA1 and A2A1 - # using cross multiplication to avoid division by zero - return cls.is_x((p.y - a1.y) * (a2.x - a1.x), (a2.y - a1.y) * (p.x - a1.x)) - + def are_edges_parallel(cls, edge1, edge2): edge1_dir = edge1[1] - edge1[0] edge2_dir = edge2[1] - edge2[0] + return cls.is_x(edge1_dir.cross(edge2_dir).length_squared, 0) - if cls.is_x(edge1_dir.cross(edge2_dir).length_squared, 0): # check they are parallel - if is_point_on_line(edge1[0], edge2) or is_point_on_line(edge1[1], edge2): - return True - return False + @classmethod + def are_edges_collinear(cls, edge1, edge2): + if not cls.are_edges_parallel(edge1, edge2): + return False + return cls.are_edges_parallel((edge2[0], edge1[0]), edge2) @classmethod def closest_points(cls, edge1, edge2) -> bool: diff --git a/src/blenderbim/test/tool/test_cad.py b/src/blenderbim/test/tool/test_cad.py new file mode 100644 index 0000000000..1aaa098fde --- /dev/null +++ b/src/blenderbim/test/tool/test_cad.py @@ -0,0 +1,58 @@ +# BlenderBIM Add-on - OpenBIM Blender Add-on +# Copyright (C) 2023 Dion Moult , @Andrej730 +# +# This file is part of BlenderBIM Add-on. +# +# BlenderBIM Add-on is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# BlenderBIM Add-on is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with BlenderBIM Add-on. If not, see . + +from test.bim.bootstrap import NewFile +from blenderbim.tool.cad import Cad as subject +from mathutils import Vector + +V = lambda *x: Vector([float(i) for i in x]) + + +class TestAreEdgesCollinear(NewFile): + def test_run(self): + # fmt: off + # Parallel edges but not collinear (different z-coordinates) + assert not subject.are_edges_collinear( + (V(-1,0,-1), V(1,0,-1)), + (V(-1,0,1), V(1,0,1)) + ) + + # One edge is just a point and the other is a line segment. + assert not subject.are_edges_collinear( + (V(1,-1,0), V(1,-1,0)), + (V(-1,1,0), V(1,1,0)) + ) + + # Both edges are collinear and overlap. + assert subject.are_edges_collinear( + (V(0,0,0), V(2,2,2)), + (V(1,1,1), V(3,3,3)) + ) + + # Both edges are collinear but don't overlap. + assert subject.are_edges_collinear( + (V(0,0,0), V(1,1,1)), + (V(2,2,2), V(3,3,3)) + ) + + # Edges are not parallel and not collinear. + assert not subject.are_edges_collinear( + (V(0,0,0), V(1,1,1)), + (V(0,1,0), V(1,0,1)) + ) + # fmt: on diff --git a/src/ifcopenshell-python/ifcopenshell/util/shape_builder.py b/src/ifcopenshell-python/ifcopenshell/util/shape_builder.py index 6cf52d0581..64a93cf211 100644 --- a/src/ifcopenshell-python/ifcopenshell/util/shape_builder.py +++ b/src/ifcopenshell-python/ifcopenshell/util/shape_builder.py @@ -19,12 +19,17 @@ import collections import ifcopenshell import ifcopenshell.api -from math import cos, sin, pi, tan, radians +from math import cos, sin, pi, tan, radians, degrees, atan, sqrt from mathutils import Vector, Matrix from itertools import chain V = lambda *x: Vector([float(i) for i in x]) sign = lambda x: x and (1, -1)[x < 0] +PRECISION = 1.0e-5 +is_x = lambda value, x: (x + PRECISION) > value > (x - PRECISION) +round_to_precision = lambda x, si_conversion: round(x * si_conversion, 5) / si_conversion +round_vector_to_precision = lambda v, si_conversion: Vector([round_to_precision(i, si_conversion) for i in v]) + # Note: using ShapeBuilder try not to reuse IFC elements in the process # otherwise you might run into situation where builder.mirror or other operation @@ -840,7 +845,9 @@ class ShapeBuilder: return face_set - def mep_transition_shape(self, start_segment, end_segment, start_length, end_length, angle=30.0): + def mep_transition_shape( + self, start_segment, end_segment, start_length, end_length, angle=30.0, profile_offset=None + ): """ returns tuple of Model/Body/MODEL_VIEW IfcRepresentation and transition shape data """ @@ -853,68 +860,285 @@ class ShapeBuilder: if material and material.is_a("IfcMaterialProfileSet") and len(material.MaterialProfiles) == 1: return material.MaterialProfiles[0].Profile + def get_circle_points(radius, segments=16): + """starting from (R,0), going counter-clockwise""" + angle_d = 2 * pi / segments + verts = [] + for i in range(segments): + angle = angle_d * i + verts.append(V(cos(angle), sin(angle), 0) * radius) + return verts + + def get_rectangle_points(dim): + """Starting from (+X/2, +Y/2) going counter-clockwise""" + dim = dim / 2 + points = [ + dim * V(1, 1, 0), + dim * V(-1, 1, 0), + dim * V(-1, -1, 0), + dim * V(1, -1, 0), + ] + return points + + # TODO: support more profiles + def get_dim(profile, depth): + if profile.is_a("IfcRectangleProfileDef"): + return V(profile.XDim / 2, profile.YDim / 2, depth) + elif profile.is_a("IfcCircleProfileDef"): + return V(profile.Radius, profile.Radius, depth) + return None + + def get_profile_faceset(points, length, offset=None): + # prevent mutating arguments, deepcopy doesn't work + start_points = [p.copy() if not offset else (p + offset) for p in points] + end_points = [p.copy() for p in start_points] + for p in end_points: + p.z += length + + points = start_points + end_points + faces = [] + n_verts = len(start_points) + last_vert_i = n_verts - 1 + for i in range(last_vert_i): + face = (i, i + 1, n_verts + i + 1, n_verts + i) + faces.append(face) + faces.append((last_vert_i, 0, n_verts + 0, n_verts + last_vert_i)) # close the loop + + # if there is offset we put a cap at the end + # otherwise at the start + if offset: + faces.append(tuple(range(n_verts, n_verts * 2))) + else: + faces.append(tuple(reversed(range(n_verts)))) + + face_set = self.polygonal_face_set(points, faces) + return face_set + start_profile = get_profile(start_segment) end_profile = get_profile(end_segment) - # TODO: support more profiles - if not start_profile.is_a("IfcRectangleProfileDef") or not end_profile.is_a("IfcRectangleProfileDef"): - # Non rectangular profiles are not yet supported + start_half_dim = get_dim(start_profile, start_length) + end_half_dim = get_dim(end_profile, end_length) + + # if profile types are not supported + if not start_half_dim or not end_half_dim: return None, None - start_half_dim = V(start_profile.XDim / 2, start_profile.YDim / 2, start_length) - end_half_dim = V(end_profile.XDim / 2, end_profile.YDim / 2, end_length) - transition_items = [] - end_extrusion_offset = V(0, 0, start_length) + start_offset = V(0, 0, start_length) + end_extrusion_offset = start_offset.copy() + + si_conversion = ifcopenshell.util.unit.calculate_unit_scale(self.file) + + # TODO: support offseted profiles + def get_transition_length(start_half_dim, end_half_dim, angle, profile_offset=None): + # NOTE: transition_length == 0 when profiles have the same dimensions + + # holy grail of the transition length: - def get_transition_legth(start_half_dim, end_half_dim, angle): diff = start_half_dim.xy - end_half_dim.xy diff = Vector([abs(i) for i in diff]) - c = diff.x * tan(radians(90 - angle / 2)) - a = diff.y - b = (c**2 - a**2) ** 0.5 - return b - transition_length = get_transition_legth(start_half_dim, end_half_dim, angle) + def calculate_transition(diff, profile_offset, end_profile=False, angle=None, length=None): + """will return transition length based on the profile dimension differences and offset. + + If `length` is provided will return transition angle""" + + # offsets tend to have bunch of float point garbage + # that can result in errors when we're calculating value for square root below + offset = V(0, 0) if profile_offset is None else round_vector_to_precision(profile_offset, si_conversion) + if end_profile: + diff, offset = diff.yx, offset.yx + + a = diff.x + offset.x + b = diff.x - offset.x + if length is None: + if diff.x == 0: + return 0 + + t = tan(radians(angle)) + l1 = (a + b + sqrt(a**2 + 4 * a * b * t**2 + 2 * a * b + b**2)) / (2 * t) + length = sqrt(l1**2 - offset.y**2) + + # TODO: remove after debug, move somewhere to tests? + if True: + A = (end_profile if end_profile else start_half_dim) * V(1, 0, 0) + end_profile_offset = offset.to_3d() + V(0, 0, length) + D = (start_half_dim if end_profile else end_half_dim) * V(1, 0, 0) + B, C = -A, -D + C += end_profile_offset + D += end_profile_offset + tested_angle = degrees((A - D).angle(B - C)) + print(f"length = {length}, requested angle = {angle}, tested angle = {tested_angle}") + return length + + elif angle is None: + # TODO: need to handle angle differently for that case + # it occurs when diff == 0 + if length == 0: + return 0 + + l1 = sqrt(length**2 + offset.y**2) + t = -l1 * (a + b) / (a * b - l1**2) + angle = atan(t) + return angle + + transition_lengths = [ + calculate_transition(diff, profile_offset, angle=angle), + calculate_transition(diff, profile_offset, angle=angle, end_profile=True), + ] + + other_side_angles = [ + calculate_transition(diff, profile_offset, length=transition_lengths[0]), + calculate_transition(diff, profile_offset, length=transition_lengths[1], end_profile=True), + ] + + # NOTE: debug values + print(f"offset = {profile_offset}") + print(f"diff = {diff}") + print(f"lengths = {transition_lengths}") + print(f"other angles = {other_side_angles}") + print(f"measurable angles = {[(180 - deg)/2 for deg in other_side_angles]}") + + # need to make sure that the worst angle (maximum angle) + # for this transition angle is `angle` + for transition_length, other_side_angle in zip(transition_lengths, other_side_angles): + if other_side_angle < angle or is_x(other_side_angle, angle): + print(f"final length = {transition_length}") # TODO: remove after debug + return transition_length + + transition_length = get_transition_length(start_half_dim, end_half_dim, angle, profile_offset) + if transition_length is None: + return None, None + faces = [] - if transition_length != 0: - end_extrusion_offset.z += transition_length + end_extrusion_offset.z += transition_length + if profile_offset: + end_extrusion_offset.xy += profile_offset + + if start_profile.is_a("IfcRectangleProfileDef") and end_profile.is_a("IfcRectangleProfileDef"): + # no transitions for exactly the same profiles + if transition_length == 0: + return None, None faces += [(3, 4, 7, 0), (11, 8, 15, 12), (3, 11, 12, 4), (7, 15, 8, 0)] - # NOTE: clockwise order for correct face orientation - faces += [ - # start extrusion - (0, 1, 2, 3), - (8, 11, 10, 9), - (0, 8, 9, 1), - (1, 9, 10, 2), - (2, 10, 11, 3), - # end extrusion - (4, 5, 6, 7), - (12, 15, 14, 13), - (4, 12, 13, 5), - (5, 13, 14, 6), - (6, 14, 15, 7), - ] - points = [ - start_half_dim * V(-1, -1, 1), - start_half_dim * V(-1, -1, 0), - start_half_dim * V(1, -1, 0), - start_half_dim * V(1, -1, 1), - end_half_dim * V(1, -1, 0) + end_extrusion_offset, - end_half_dim * V(1, -1, 1) + end_extrusion_offset, - end_half_dim * V(-1, -1, 1) + end_extrusion_offset, - end_half_dim * V(-1, -1, 0) + end_extrusion_offset, - start_half_dim * V(-1, 1, 1), - start_half_dim * V(-1, 1, 0), - start_half_dim * V(1, 1, 0), - start_half_dim * V(1, 1, 1), - end_half_dim * V(1, 1, 0) + end_extrusion_offset, - end_half_dim * V(1, 1, 1) + end_extrusion_offset, - end_half_dim * V(-1, 1, 1) + end_extrusion_offset, - end_half_dim * V(-1, 1, 0) + end_extrusion_offset, - ] + # NOTE: clockwise order for correct face orientation + faces += [ + # start extrusion + (0, 1, 2, 3), + (8, 11, 10, 9), + (0, 8, 9, 1), + (1, 9, 10, 2), + (2, 10, 11, 3), + # end extrusion + (4, 5, 6, 7), + (12, 15, 14, 13), + (4, 12, 13, 5), + (5, 13, 14, 6), + (6, 14, 15, 7), + ] + points = [ + start_half_dim * V(-1, -1, 1), + start_half_dim * V(-1, -1, 0), + start_half_dim * V(1, -1, 0), + start_half_dim * V(1, -1, 1), + end_half_dim * V(1, -1, 0) + end_extrusion_offset, + end_half_dim * V(1, -1, 1) + end_extrusion_offset, + end_half_dim * V(-1, -1, 1) + end_extrusion_offset, + end_half_dim * V(-1, -1, 0) + end_extrusion_offset, + start_half_dim * V(-1, 1, 1), + start_half_dim * V(-1, 1, 0), + start_half_dim * V(1, 1, 0), + start_half_dim * V(1, 1, 1), + end_half_dim * V(1, 1, 0) + end_extrusion_offset, + end_half_dim * V(1, 1, 1) + end_extrusion_offset, + end_half_dim * V(-1, 1, 1) + end_extrusion_offset, + end_half_dim * V(-1, 1, 0) + end_extrusion_offset, + ] + elif start_profile.is_a("IfcCircleProfileDef") and end_profile.is_a("IfcCircleProfileDef"): + # no transitions for exactly the same profiles + if transition_length == 0: + return None, None + + n_segments = 16 + first_profile_points = get_circle_points(start_profile.Radius, n_segments) + second_profile_points = get_circle_points(end_profile.Radius, n_segments) + + faces = [] + for i in range(n_segments): + # For wrapping around the circle + next_i = (i + 1) % n_segments + face = [i, next_i, next_i + n_segments, i + n_segments] + faces.append(face) + + transition_items.append(get_profile_faceset(first_profile_points, start_length)) + transition_items.append(get_profile_faceset(second_profile_points, end_length, end_extrusion_offset)) + + first_profile_points = [p + start_offset for p in first_profile_points] + second_profile_points = [p + end_extrusion_offset for p in second_profile_points] + + points = first_profile_points + second_profile_points + + else: # one is circular, another one is rectangular + # support transition from rectangle to circle of the same dimensions + if transition_length == 0: + transition_length = (start_length + end_length) / 2 + end_extrusion_offset.z += transition_length + + starting_with_circle = start_profile.is_a("IfcCircleProfileDef") + if starting_with_circle: + circle_profile, rect_profile = start_profile, end_profile + else: + circle_profile, rect_profile = end_profile, start_profile + + circle_points = get_circle_points(circle_profile.Radius) + rect_points = get_rectangle_points(V(rect_profile.XDim, rect_profile.YDim, 0)) + + if starting_with_circle: + start_points, end_points = circle_points, rect_points + else: + start_points, end_points = rect_points, circle_points + + transition_items.append(get_profile_faceset(start_points, start_length)) + transition_items.append(get_profile_faceset(end_points, end_length, end_extrusion_offset)) + + # offset verts + if starting_with_circle: + circle_points = [p + start_offset for p in circle_points] + rect_points = [p + end_extrusion_offset for p in rect_points] + else: + rect_points = [p + start_offset for p in rect_points] + circle_points = [p + end_extrusion_offset for p in circle_points] + + # circle verts are 0-15, rect verts are 16-19 + points = circle_points + rect_points + transition_faces = [ + (0, 19, 16), # base + (0, 16, 1), + (1, 16, 2), + (2, 16, 3), + (3, 16, 4), + (4, 16, 17), # base + (4, 17, 5), + (5, 17, 6), + (6, 17, 7), + (7, 17, 8), + (8, 17, 18), # base + (8, 18, 9), + (9, 18, 10), + (10, 18, 11), + (11, 18, 12), + (12, 18, 19), # base + (12, 19, 13), + (13, 19, 14), + (14, 19, 15), + (15, 19, 0), + ] + # revert them in case it's starting with circle profile to keep the face orientation + if starting_with_circle: + transition_faces = [f[::-1] for f in transition_faces] + faces += transition_faces face_set = self.polygonal_face_set(points, faces) transition_items.append(face_set)