bim.mep_add_bend

Added operator to create bend fittings between objects, it's still work in progress and now it only supports rectangular profiles.

Small demo - https://imgur.com/a/vM7hnJa

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

@@ -179,6 +179,7 @@ classes = (
     roof.SetGableRoofEdgeAngle,
     mep.MEPAddObstruction,
     mep.MEPAddTransition,
+    mep.MEPAddBend,
 )
 
 addon_keymaps = []

src/blenderbim/blenderbim/bim/module/model/mep.py

@@ -35,7 +35,7 @@ import blenderbim.core.type
 import blenderbim.core.root
 import blenderbim.core.geometry
 import blenderbim.tool as tool
-from math import pi, degrees, radians
+from math import pi, degrees, radians, sin, cos, asin
 from copy import copy
 from mathutils import Vector, Matrix
 from ifcopenshell.util.shape_builder import ShapeBuilder
@@ -201,6 +201,7 @@ class FitFlowSegments(bpy.types.Operator, tool.Ifc.Operator):
                 is_on_axis2 = tool.Cad.is_point_on_edge(intersect2, axis2)
                 if not is_on_axis1 and not is_on_axis2:
                     fitting_type = "BEND"
+                    bpy.ops.bim.mep_add_bend()
                 elif is_on_axis1 and is_on_axis2:
                     fitting_type = "CROSS"
                 else:
@@ -644,32 +645,40 @@ class MEPAddTransition(bpy.types.Operator, tool.Ifc.Operator):
             return {"CANCELLED"}
 
         # TODO: support different profiles rotation by local Z
-        rotation_difference_z = start_object.matrix_world.to_quaternion().rotation_difference(end_object.matrix_world.to_quaternion()).to_euler().z
+        # check rotation difference
+        end_object_rotation = end_object.matrix_world.to_quaternion()
+        rotation_difference_z = (
+            start_object.matrix_world.to_quaternion().rotation_difference(end_object_rotation).to_euler().z
+        )
+
         def is_multiple_of_pi(value):
             n = round(value / pi)
             return tool.Cad.is_x(abs(value - n * pi), 0)
-        
+
         if not is_multiple_of_pi(rotation_difference_z):
-            self.report({"ERROR"}, f"There is some rotation difference between profiles by local Z axis: {round(degrees(rotation_difference_z))} deg, this kind of transition is not yet supported.")
+            self.report(
+                {"ERROR"},
+                "There is some rotation difference between profiles by local Z axis: "
+                f"{round(degrees(rotation_difference_z))} deg, this kind of transition is not yet supported.",
+            )
             return {"CANCELLED"}
 
+        # setup start / end points
         start_segment_data = MEPGenerator().get_segment_data(start_element)
         end_segment_data = MEPGenerator().get_segment_data(end_element)
-        end_port = end_segment_data["start_port"]
-        start_port = start_segment_data["end_port"]
-
         points_ports_map = {
             start_segment_data["start_point"]: start_segment_data["start_port"],
             start_segment_data["end_point"]: start_segment_data["end_port"],
             end_segment_data["start_point"]: end_segment_data["start_port"],
             end_segment_data["end_point"]: end_segment_data["end_port"],
         }
-
         # transition points
-        start_point, end_point = tool.Cad.closest_points(
+        (start_point, end_point), (first_segment_start, second_segment_end) = tool.Cad.closest_points(
             (start_segment_data["start_point"], start_segment_data["end_point"]),
             (end_segment_data["start_point"], end_segment_data["end_point"]),
         )
+        start_port = points_ports_map[start_point]
+        end_port = points_ports_map[end_point]
 
         # figure profile offset
         base_transition_dir = keep_only_z_axis(end_point - start_point).normalized()
@@ -694,18 +703,6 @@ class MEPAddTransition(bpy.types.Operator, tool.Ifc.Operator):
             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"],
-            )
-            if p not in (start_point, end_point)
-        ]
-
         def get_segments_length():
             start_dir = (start_point - first_segment_start).normalized()
             segments_vector = second_segment_end - first_segment_start
@@ -716,8 +713,6 @@ class MEPAddTransition(bpy.types.Operator, tool.Ifc.Operator):
         # 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)
@@ -752,6 +747,8 @@ class MEPAddTransition(bpy.types.Operator, tool.Ifc.Operator):
         # adjust the segments
         end_object_rotation = end_object.matrix_world.to_quaternion()
         end_object_z_basis = end_object_rotation.to_matrix().col[2]  # z basis vector
+
+        # TODO: do it beforehand, as with bends
         if tool.Cad.is_x(start_object_z_basis.dot(transition_dir), 1):
             start_connection = "ATEND"
         else:
@@ -768,13 +765,11 @@ class MEPAddTransition(bpy.types.Operator, tool.Ifc.Operator):
             [start_element, end_element], [start_port, end_port], "TRANSITION"
         )
         transition_type = fitting_data["fitting_type"] if fitting_data else None
+        start_port_match = fitting_data["start_port_match"] if fitting_data else True
         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
-
-        # create new fitting type if nothing is compatible
-        if not transition_type:
+        else:  # create new fitting type if nothing is compatible
             mesh = bpy.data.meshes.new("Transition")
             obj = bpy.data.objects.new("Transition", mesh)
             transition_type = blenderbim.core.root.assign_class(
@@ -825,3 +820,316 @@ class MEPAddTransition(bpy.types.Operator, tool.Ifc.Operator):
         tool.Ifc.run("system.connect_port", port1=ports[1], port2=end_port, direction="NOTDEFINED")
 
         return {"FINISHED"}
+
+
+class MEPAddBend(bpy.types.Operator, tool.Ifc.Operator):
+    bl_idname = "bim.mep_add_bend"
+    bl_label = "Add Bend"
+    bl_description = "Adds a bend between two MEP elements. Elements are either provided by ID or selected in Blender"
+    bl_options = {"REGISTER", "UNDO"}
+    start_length: bpy.props.FloatProperty(
+        name="Start Length", description="Bend start length in SI units", default=0.1, subtype="DISTANCE"
+    )
+    end_length: bpy.props.FloatProperty(
+        name="End Length", description="Bend end length in SI units", default=0.1, subtype="DISTANCE"
+    )
+    start_segment_id: bpy.props.IntProperty(name="Start Segment Element ID", default=0)
+    end_segment_id: bpy.props.IntProperty(name="End Segment Element ID", default=0)
+    radius: bpy.props.FloatProperty(
+        "Bend Inner Radius", description="Bend inner radius in SI units", default=0.2, subtype="DISTANCE"
+    )
+
+    def _execute(self, context):
+        start_element, end_element = None, None
+        ifc_file = tool.Ifc.get()
+        si_conversion = ifcopenshell.util.unit.calculate_unit_scale(ifc_file)
+
+        self.start_length, self.end_length = 0, 0
+
+        if not (self.start_length == 0 and self.end_length == 0):
+            self.report({"ERROR"}, f"Only zero lengths are now supported.")
+            return {"CANCELLED"}
+
+        if self.start_segment_id and self.end_segment_id:
+            start_element = ifc_file.by_id(self.start_segment_id)
+            end_element = ifc_file.by_id(self.end_segment_id)
+            start_object = tool.Ifc.get_object(start_element)
+            end_object = tool.Ifc.get_object(end_element)
+
+        elif len(context.selected_objects) == 2:
+            start_object = context.active_object
+            end_object = next(o for o in context.selected_objects if o != context.active_object)
+            start_element = tool.Ifc.get_entity(start_object)
+            end_element = tool.Ifc.get_entity(end_object)
+            if not start_element or not end_element:
+                self.report({"ERROR"}, f"Two IFC elements should be selected for the bend.")
+                return {"CANCELLED"}
+
+        else:
+            self.report({"ERROR"}, f"Two IFC elements should be provided for the bend.")
+            return {"CANCELLED"}
+
+        # check rotation difference
+        def rotation_difference_check():
+            end_object_rotation = end_object.matrix_world.to_quaternion()
+            rotation_difference = start_object.matrix_world.to_quaternion().rotation_difference(end_object_rotation).to_euler()
+
+            def is_multiple_of_pi(value):
+                n = round(value / pi)
+                return tool.Cad.is_x(abs(value - n * pi), 0)
+
+            if not is_multiple_of_pi(rotation_difference.z):
+                error_msg = (
+                    "There is some rotation difference between profiles by local Z axis: "
+                    f"{round(degrees(rotation_difference.z))} deg, adding a bend is not possible."
+                )
+                return error_msg
+
+        if error_msg := rotation_difference_check():
+            self.report({"ERROR"}, error_msg)
+            return {"CANCELLED"}
+
+        # check segments types
+        def types_check():
+            start_type = ifcopenshell.util.element.get_type(start_element)
+            end_type = ifcopenshell.util.element.get_type(end_element)
+            if not start_type or not end_type:
+                return False
+            return start_type == end_type
+
+        if not types_check():
+            self.report(
+                {"ERROR"},
+                "Segments types do not match "
+                "or one of the segments doesn't have type which is required for a bend.",
+            )
+            return {"CANCELLED"}
+
+        # TODO: support circular profiles
+        profile = tool.Model.get_flow_segment_profile(start_element)
+        if not profile.is_a("IfcRectangleProfileDef"):
+            self.report(
+                {
+                    "ERROR",
+                    "For now Only IfcRectangleProfileDef profiles supported for a bend, "
+                    f"the segments are {profile.is_a()}",
+                }
+            )
+            return {"CANCELLED"}
+
+        def get_dim(profile):
+            if profile.is_a("IfcRectangleProfileDef"):
+                return V(profile.XDim / 2, profile.YDim / 2)
+            elif profile.is_a("IfcCircleProfileDef"):
+                return V(profile.Radius, profile.Radius)
+            return None
+
+        # setup start / end points
+        start_object_rotation = start_object.matrix_world.to_quaternion().to_matrix()
+        start_segment_data = MEPGenerator().get_segment_data(start_element)
+        end_segment_data = MEPGenerator().get_segment_data(end_element)
+        points_ports_map = {
+            start_segment_data["start_point"]: start_segment_data["start_port"],
+            start_segment_data["end_point"]: start_segment_data["end_port"],
+            end_segment_data["start_point"]: end_segment_data["start_port"],
+            end_segment_data["end_point"]: end_segment_data["end_port"],
+        }
+        (start_point, end_point), (first_segment_start, second_segment_end) = tool.Cad.closest_points(
+            (start_segment_data["start_point"], start_segment_data["end_point"]),
+            (end_segment_data["start_point"], end_segment_data["end_point"]),
+        )
+        start_port = points_ports_map[start_point]
+        end_port = points_ports_map[end_point]
+        start_point_on_origin = start_point == start_segment_data["start_point"]
+        start_connection = "ATSTART" if start_point_on_origin else "ATEND"
+        start_segment_sign = -1 if start_point_on_origin else 1
+
+        end_point_on_origin = end_point == end_segment_data["start_point"]
+        end_connection = "ATSTART" if end_point_on_origin else "ATEND"
+        end_segment_sign = -1 if end_point_on_origin else 1
+
+        profile_dim = get_dim(profile) * si_conversion
+
+        # TODO: profile offset may need to be flipped (check transition code)
+        to_start_object_space = start_object_rotation.inverted()
+        profile_offset = (to_start_object_space @ end_point) - (to_start_object_space @ start_point)
+
+        def check_for_double_bends():
+            # The theory is To avoid double bends, the profile offset should occur along only two axes:
+            # 1) The local Z-axis of the start segment
+            # 2) One of the lateral axes (either X or Y)
+            #
+            # Double bend required when:
+            # - there are 2 or 0 lateral axes involved
+            # - offset appear by the non-lateral axis
+            #
+            # NOTE: some double bends are only possible for square profiles:
+            # https://i.imgur.com/ZhdGbEp.png
+
+            z_axis_end_object = end_object.matrix_world.col[2].normalized().to_3d()
+            z_axis_end_object_local = to_start_object_space @ z_axis_end_object
+            lateral_axes = [i for i in range(2) if not tool.Cad.is_x(z_axis_end_object_local[i], 0)]
+
+            if len(lateral_axes) != 1:
+                return (
+                    None,
+                    f"For now only one lateral axis is supported for a bend (double bends not supported). Found lateral axes: {len(lateral_axes)}.",
+                )
+
+            non_lateral_axis = 0 if lateral_axes[0] == 1 else 1
+            non_lateral_axis_offset = profile_offset[non_lateral_axis]
+            if not tool.Cad.is_x(non_lateral_axis_offset, 0):
+                return (
+                    None,
+                    "For now offset by non-lateral axis is not supported for a bend (double bends not supported).\n"
+                    f"Detected an offset of {round(non_lateral_axis_offset, 5)} along the local axis {'XY'[non_lateral_axis]} when lateral axis is {'XY'[lateral_axes[0]]}.",
+                )
+
+            return lateral_axes[0], None
+
+        lateral_axis, error_msg = check_for_double_bends()
+        if error_msg:
+            self.report({"ERROR"}, error_msg)
+            return {"CANCELLED"}
+
+        O = V(0, 0, 0)
+        get_z_basis = lambda o: o.matrix_world.col[2].normalized().to_3d()
+        angle = tool.Cad.angle_edges((get_z_basis(start_object), O), (get_z_basis(end_object), O))
+
+        lateral_sign = tool.Cad.sign(profile_offset[lateral_axis])
+        radial_offset = V(0, 0, 0)
+        ref_point_radius = self.radius + profile_dim[lateral_axis]
+        radial_offset[lateral_axis] = ref_point_radius * (1 - cos(angle)) * lateral_sign
+        radial_offset.z = ref_point_radius * sin(angle)
+
+        def get_segments_extend():
+            end_segment_z_local = to_start_object_space @ get_z_basis(end_object)
+            segments_intersection = tool.Cad.intersect_edges(
+                (V(0, 0, 1), V(0, 0, 0)), (profile_offset + end_segment_z_local, profile_offset)
+            )[0]
+
+            curent_start_offset = segments_intersection.length
+            required_start_offset = abs(radial_offset.z)
+            current_end_offset = (segments_intersection - profile_offset).length
+            required_end_offset = abs(radial_offset[lateral_axis])
+
+            start_extend = curent_start_offset - required_start_offset
+            end_extend = current_end_offset - required_end_offset
+
+            return start_extend, end_extend
+
+        def check_new_segment_length(start_point, end_point, extend_point):
+            """Check if segment is placed too near to the bend point.
+
+            The idea is that we can either extend segment toward the bend
+            but we can shrink it only until it's start.
+
+            If the segment is too near it will return offset to fix the problem,
+            otherwise returns `None`.
+
+            """
+            base_edge = end_point - start_point
+            new_edge = extend_point - start_point
+            projection = new_edge.dot(base_edge.normalized())
+            if projection < 0 or tool.Cad.is_x(projection, 0):
+                return projection
+            return None
+
+        # adjust segments to fit the radius and angle
+        start_segment_extend, end_segment_extend = get_segments_extend()
+
+        start_segment_extend_point = start_point + start_segment_sign * start_segment_extend * get_z_basis(start_object)
+        projection = check_new_segment_length(first_segment_start, start_point, start_segment_extend_point)
+        if projection is not None:
+            self.report(
+                {"ERROR"},
+                f"Start segment starts too near to the bend, need to offset it atleast by {round(projection, 3)} m.",
+            )
+            return {"ERROR"}
+
+        end_segment_extend_point = end_point + end_segment_sign * end_segment_extend * get_z_basis(end_object)
+        projection = check_new_segment_length(second_segment_end, end_point, end_segment_extend_point)
+        if projection is not None:
+            self.report(
+                {"ERROR"},
+                f"End segment starts too near to the bend, need to offset it atleast by {round(projection, 3)} m.",
+            )
+            return {"ERROR"}
+
+        DumbProfileJoiner().join_E(start_object, start_segment_extend_point, start_connection)
+        DumbProfileJoiner().join_E(end_object, end_segment_extend_point, end_connection)
+
+        context.view_layer.update()  # update matrices
+
+        builder = ShapeBuilder(ifc_file)
+        rep, bend_data = builder.mep_bend_shape(
+            start_element,
+            self.start_length / si_conversion,
+            self.end_length / si_conversion,
+            angle,
+            self.radius / si_conversion,
+            profile_offset / si_conversion,
+        )
+
+        bpy.ops.bim.create_shape_from_step_id(step_id=rep.id(), should_include_curves=True)
+
+        # find the compatible fitting type
+        fitting_data = MEPGenerator().get_compatible_fitting_type(
+            [start_element, end_element], [start_port, end_port], "BEND"
+        )
+        bend_type = fitting_data["fitting_type"] if fitting_data else None
+        start_port_match = fitting_data["start_port_match"] if fitting_data else True
+        if bend_type:
+            # TODO: handle the case without creating a representation in the first place?
+            ifcopenshell.api.run("geometry.remove_representation", ifc_file, representation=rep)
+        else:  # create new fitting type if nothing is compatible
+            mesh = bpy.data.meshes.new("Bend")
+            obj = bpy.data.objects.new("Bend", mesh)
+            bend_type = blenderbim.core.root.assign_class(
+                tool.Ifc,
+                tool.Collector,
+                tool.Root,
+                obj=obj,
+                ifc_class=MEPGenerator().get_mep_element_class_name(start_element, "FittingType"),
+                predefined_type="BEND",
+                should_add_representation=False,
+            )
+            body = ifcopenshell.util.representation.get_context(ifc_file, "Model", "Body", "MODEL_VIEW")
+            tool.Model.replace_object_ifc_representation(body, obj, rep)
+            pset = ifcopenshell.api.run("pset.add_pset", tool.Ifc.get(), product=bend_type, name="BBIM_Fitting")
+            ifcopenshell.api.run(
+                "pset.edit_pset",
+                tool.Ifc.get(),
+                pset=pset,
+                properties={"Data": json.dumps(bend_data, default=list)},
+            )
+
+        # NOTE: at this point we loose current blender objects selection
+        # create transition element
+        bpy.ops.bim.add_constr_type_instance(relating_type_id=bend_type.id())
+        fitting_obj = bpy.context.active_object
+
+        # adjust fitting object rotation and location
+        # required since we'll base our `fitting_obj_dir` on this
+        fitting_obj.matrix_world = start_object.matrix_world
+        context.view_layer.update()
+
+        # depending on fitting direction we may need to flip it or attach it's origin to end segment
+        # direction can be different depending on:
+        # - order of the current segments
+        # - order of the segments that were used with the same fitting type before
+        direction_match = tool.Cad.are_vectors_equal(get_z_basis(start_object), get_z_basis(fitting_obj))
+        # if there are no mismatches or everything matches up we don't need to flip the transition
+        if start_port_match != direction_match:
+            fitting_obj.matrix_world = start_object.matrix_world @ Matrix.Rotation(radians(180), 4, "X")
+        fitting_obj.location = start_segment_extend_point if start_port_match else end_segment_extend_point
+
+        # add ports and connect them
+        ports = tool.System.add_ports(fitting_obj, offset_end_port=start_object_rotation @ (radial_offset * V(1, 1, 0)))
+        if not start_port_match:
+            start_port, end_port = end_port, start_port
+        tool.Ifc.run("system.connect_port", port1=ports[0], port2=start_port, direction="NOTDEFINED")
+        tool.Ifc.run("system.connect_port", port1=ports[1], port2=end_port, direction="NOTDEFINED")
+
+        self.report({"INFO"}, f"Success!.. kind of. The angle was {round(bend_data['angle'])}")
+        return {"FINISHED"}

src/blenderbim/blenderbim/tool/cad.py

@@ -35,6 +35,7 @@ import math
 import bmesh
 import mathutils.geometry
 from mathutils import Vector, Matrix, geometry
+import itertools
 
 
 VTX_PRECISION = 1.0e-5
@@ -74,6 +75,8 @@ class Cad:
         """
         > takes 2 edges, each as a tuple of two vectors
         < returns the potentially signed angle as degrees or radians
+
+        NOTE: `signed` expects both edges to be 2D (just as `Vector.angle_signed`)
         """
         if signed:
             a = (edge1[1] - edge1[0]).angle_signed(edge2[1] - edge2[0])
@@ -248,23 +251,33 @@ class Cad:
         return cls.are_edges_parallel((edge2[0], edge1[0]), edge2)
 
     @classmethod
-    def closest_points(cls, edge1, edge2) -> bool:
+    def closest_points(cls, edge1, edge2):
+        """
+        closest end points between `edge1` and `edge2`
+
+        ensures returned vectors are the exact objects
+        that were passed to the method with `edge1` and `edge2`
+
+        < returns two tuples - two closest points and two other points
+
+        first point of each tuple belongs to `edge1` and second to `edge2`
         """
 
-        closest end points between `edge1` and `edge2` assuming `edge1` and `edge2` are collinear.
+        distance_squared = None
+        closest_points = None
+        for p1 in edge1:
+            for p2 in edge2:
+                cur_line = p2 - p1
+                cur_distance_squared = cur_line.dot(cur_line)
+                if distance_squared is None or cur_distance_squared < distance_squared:
+                    closest_points = (p1, p2)
+                    distance_squared = cur_distance_squared
 
-        < returns two points, first one belongs to `edge1` and second to `edge2`
-
-        """
-        direction = (edge1[1] - edge1[0]).normalized()
-
-        # Project points onto the line to get scalar values along the direction
-        points_values = [(p, p.dot(direction)) for p in (edge1 + edge2)]
-        sorted_points = sorted(points_values, key=lambda el: el[1])
-
-        edge1_point = next((p for p, v in sorted_points[1:3] if p in edge1), None)
-        edge2_point = next((p for p, v in sorted_points[1:3] if p in edge2), None)
-        return edge1_point, edge2_point
+        other_points = (
+            edge1[0] if closest_points[0] == edge1[1] else edge1[1],
+            edge2[0] if closest_points[1] == edge2[1] else edge2[1],
+        )
+        return closest_points, other_points
 
     @classmethod
     def find_intersecting_edges(cls, bm, pt, idx1, idx2):
@@ -489,3 +502,15 @@ class Cad:
     def is_counter_clockwise_order(cls, A, B, C):
         """whether A-B-C located in counter-clockwise order in 2d space"""
         return (C.y - A.y) * (B.x - A.x) > (B.y - A.y) * (C.x - A.x)
+
+    @classmethod
+    def sign(cls, value):
+        """
+        returns:
+          0 if cls.is_x(value, 0)) \n
+          1 if value > 0 \n
+         -1 if value < 0
+        """
+        if cls.is_x(value, 0):
+            return 0
+        return 1 if value > 0 else -1

src/blenderbim/blenderbim/tool/system.py

@@ -42,7 +42,8 @@ class System(blenderbim.core.tool.System):
             blenderbim.core.geometry.edit_object_placement(tool.Ifc, tool.Geometry, tool.Surveyor, obj=obj)
 
         mep_element = tool.Ifc.get_entity(obj)
-        length = obj.dimensions.z
+        bbox = tool.Blender.get_object_bounding_box(obj)
+        length = bbox["min_z"] if tool.Cad.is_x(bbox["max_z"], 0) else bbox["max_z"]
         ports = []
         if add_start_port:
             ports.append(add_port(mep_element, obj.matrix_world @ Matrix()))

src/blenderbim/test/tool/test_cad.py

@@ -56,3 +56,34 @@ class TestAreEdgesCollinear(NewFile):
             (V(0,1,0), V(1,0,1))
         )
         # fmt: on
+
+
+class TestClosestPoints(NewFile):
+    def test_run(self):
+        # non collinear
+        edge1 = (V(0, 0, 0), V(1, 0, 0))
+        edge2 = (V(2, 0, 1), V(2, 0, 2))
+        assert subject.closest_points(edge1, edge2)[0] == (edge1[1], edge2[0])
+
+        # check other points
+        assert subject.closest_points(edge1, edge2)[1] == (edge1[0], edge2[1])
+
+        # collinear
+        edge1 = (V(0, 0, 0), V(1, 0, 0))
+        edge2 = (V(3, 0, 0), V(2, 0, 0))
+        assert subject.closest_points(edge1, edge2)[0] == (edge1[1], edge2[1])
+
+        # parallel
+        edge1 = (V(0, 0, 0), V(1, 0, 0))
+        edge2 = (V(-5, 0, 0), V(-1, 0, 0))
+        assert subject.closest_points(edge1, edge2)[0] == (edge1[0], edge2[1])
+
+        # overlapping
+        edge1 = (V(0, 0, 0), V(3, 0, 0))
+        edge2 = (V(2, 0, 0), V(5, 0, 0))
+        assert subject.closest_points(edge1, edge2)[0] == (edge1[1], edge2[0])
+
+        # edge as a point
+        edge1 = (V(0, 0, 0), V(0, 0, 0))
+        edge2 = (V(1, 0, 1), V(2, 0, 2))
+        assert subject.closest_points(edge1, edge2)[0] == (edge1[0], edge2[0])

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

@@ -1,5 +1,5 @@
 # IfcOpenShell - IFC toolkit and geometry engine
-# Copyright (C) 2022 @Andrej730
+# Copyright (C) 2022, 2023 @Andrej730
 #
 # This file is part of IfcOpenShell.
 #
@@ -19,14 +19,21 @@
 import collections
 import ifcopenshell
 import ifcopenshell.api
-from math import cos, sin, pi, tan, radians, degrees, atan, sqrt
+from math import cos, sin, pi, tan, radians, degrees, atan, sqrt, ceil
 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)
+
+
+def is_x(value, x, si_conversion=None):
+    if si_conversion:
+        value = value * si_conversion
+    return (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])
 
@@ -101,7 +108,7 @@ class ShapeBuilder:
             if len(segment) == 3:
                 ifc_segments.append(self.file.createIfcArcIndex(segment))
 
-        # NOTE: IfcIndexPolyCurve support only consequtive segments
+        # NOTE: IfcIndexPolyCurve support only consecutive segments
         ifc_curve = self.file.createIfcIndexedPolyCurve(Points=ifc_points, Segments=ifc_segments)
         return ifc_curve
 
@@ -534,6 +541,9 @@ class ShapeBuilder:
 
         Position and position axes are in world space, extrusion vector in placement space defined by
         position_x_axis/position_y_axis/position_z_axis
+
+        NOTE: changing position also changes the resulting geometry origin.
+
         """
         # > profile_or_curve
         # > extrusion vector - as defined in coordinate system position_x_axis+position_z_axis
@@ -575,9 +585,9 @@ class ShapeBuilder:
         disk_solid = self.file.createIfcSweptDiskSolid(Directrix=path_curve, Radius=radius)
         return disk_solid
 
-    def get_representation(self, context, items, representation_type:str = None):
+    def get_representation(self, context, items, representation_type: str = None):
         """Create IFC representation for the specified context and items.
-        
+
         :param context: IfcGeometricRepresentationSubContext
         :param items: could be a list or single curve/IfcExtrudedAreaSolid
         :param representation_type: Explicitly specified RepresentationType, defaults to `None`.
@@ -615,8 +625,11 @@ class ShapeBuilder:
 
     # UTILITIES
     def extrude_by_y_kwargs(self):
-        """shortcut for `ShapeBuilder.extrude` to extrude by y axis.
-        it assumes you have 2d profile in xz plane and trying to extrude it by y axis"""
+        """Shortcut to get kwargs for `ShapeBuilder.extrude` to extrude by Y axis.
+
+        It assumes you have 2D profile in XZ plane and trying to extrude it by Y axis.
+
+        Extruding by Y using other kwargs might break ValidExtrusionDirection."""
         return {
             "position_x_axis": Vector((1, 0, 0)),
             "position_z_axis": Vector((0, -1, 0)),
@@ -856,6 +869,53 @@ class ShapeBuilder:
 
         return face_set
 
+    def extrude_face_set(
+        self, points, magnitude: float, extrusion_vector=V(0, 0, 1).freeze(), offset=None, start_cap=True, end_cap=True
+    ):
+        """
+        Method to extrude by creating face sets rather than creating IfcExtrudedAreaSolid.
+
+        Useful if your representation is already using face sets and you need to avoid using SweptSolid
+        to assure CorrectItemsForType.
+
+        :param points: list of points, assuming they form consecutive closed polyline.
+        :param magnitude: extrusion magnitude
+        :param type: float
+        :param extrusion_vector: extrusion direction, by default it's extruding by Z+ axis
+        :param type: Vector, optional
+        :param offset: offset from the points
+        :param type: Vector, optional
+        :param start_cap: if True, create start cap, by default it's True
+        :param type: bool, optional
+        :param end_cap: if True, create end cap, by default it's True
+        :param type: bool, optional
+
+        :return: IfcPolygonalFaceSet
+        """
+
+        # prevent mutating arguments, deepcopy doesn't work
+        start_points = [p.copy() if not offset else (p + offset) for p in points]
+        extrusion_offset = magnitude * extrusion_vector
+        end_points = [p + extrusion_offset for p in start_points]
+
+        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 end_cap:
+            faces.append(tuple(range(n_verts, n_verts * 2)))
+        if start_cap:
+            faces.append(tuple(reversed(range(n_verts))))
+
+        face_set = self.polygonal_face_set(points, faces)
+        return face_set
+
+    # TODO: move MEP to separate shape builder sub module
     def mep_transition_shape(
         self, start_segment, end_segment, start_length, end_length, angle=30.0, profile_offset=None
     ):
@@ -899,32 +959,6 @@ class ShapeBuilder:
                 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)
 
@@ -1004,8 +1038,10 @@ class ShapeBuilder:
                 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))
+            transition_items.append(self.extrude_face_set(first_profile_points, start_length, end_cap=False))
+            transition_items.append(
+                self.extrude_face_set(second_profile_points, end_length, end_extrusion_offset, start_cap=False)
+            )
 
             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]
@@ -1032,8 +1068,10 @@ class ShapeBuilder:
             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))
+            transition_items.append(self.extrude_face_set(start_points, start_length, end_cap=False))
+            transition_items.append(
+                self.extrude_face_set(end_points, end_length, end_extrusion_offset, start_cap=False)
+            )
 
             # offset verts
             if starting_with_circle:
@@ -1237,3 +1275,97 @@ class ShapeBuilder:
                 else:
                     angle = degrees(atan(offset.x / h))
             return angle
+
+    def mep_bend_shape(
+        self, segment, start_length: float, end_length: float, angle: float, radius: float, profile_offset: Vector
+    ):
+        """
+
+        :param segment: IfcFlowSegment for a bend.
+        Note that for a bend start and end segments types should match.
+
+        :param angle: bend angle, in radians
+        :param type: float
+        :param radius: bend radius
+        :param type: float
+        :param profile_offset: offset between start and end segments in local space of start segment
+            used mainly to determine the bend axes and their direction.
+            Values themselves are replaced by the radius.
+        :param type: Vector
+
+        :return: tuple of Model/Body/MODEL_VIEW IfcRepresentation and transition shape data
+        """
+
+        def get_profile(element):
+            material = ifcopenshell.util.element.get_material(element, should_skip_usage=True)
+            if material and material.is_a("IfcMaterialProfileSet") and len(material.MaterialProfiles) == 1:
+                return material.MaterialProfiles[0].Profile
+
+        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
+
+        # TODO: test with 0 radius
+        si_conversion = ifcopenshell.util.unit.calculate_unit_scale(self.file)
+        profile = get_profile(segment)
+        profile_dim = get_dim(profile, start_length)
+
+        rounded_offset = round_vector_to_precision(profile_offset, si_conversion)
+        lateral_axis = next(i for i in range(2) if not is_x(rounded_offset[i], 0))
+        non_lateral_axis = 1 if lateral_axis == 0 else 0
+        lateral_sign = sign(profile_offset[lateral_axis])
+        z_sign = sign(profile_offset.z)
+
+        rep_items = []
+
+        # bend circle center
+        O = V(0, 0, 0)
+        O[lateral_axis] = (radius + profile_dim[lateral_axis]) * lateral_sign
+
+        theta = angle
+
+        def get_circle_extrusion():
+            # get as much segment_length segments as possible
+            segment_length = pi / 20
+            num_segments = ceil(theta / segment_length)
+            theta_segments = [i * segment_length for i in range(num_segments)]
+            if not is_x(theta_segments[-1], theta):
+                theta_segments.append(theta)
+
+            inner_points, outer_points = [], []
+            r = radius
+
+            for cur_theta in theta_segments:
+                cur_theta -= pi / 2
+                inner = V(0, 0, 0)
+                # fmt: off
+                inner.z             = z_sign       * cos(cur_theta) * r
+                inner[lateral_axis] = lateral_sign * sin(cur_theta) * r
+                inner_points.append(inner)
+
+                outer = V(0, 0, 0)
+                outer.z             = z_sign       * cos(cur_theta) * (r + 2 * profile_dim[lateral_axis])
+                outer[lateral_axis] = lateral_sign * sin(cur_theta) * (r + 2 * profile_dim[lateral_axis])
+                outer_points.append(outer)
+                # fmt: on
+
+            points = inner_points + outer_points[::-1]
+            points = [p + O for p in points]
+            offset = V(0, 0, 0)
+            offset[non_lateral_axis] = -profile_dim[non_lateral_axis]
+            extrusion_vector = V(0, 0, 0)
+            extrusion_vector[non_lateral_axis] = 1
+            extrusion = self.extrude_face_set(
+                points, magnitude=profile_dim[non_lateral_axis] * 2, offset=offset, extrusion_vector=extrusion_vector
+            )
+            return extrusion
+
+        rep_items.append(get_circle_extrusion())
+        body = ifcopenshell.util.representation.get_context(self.file, "Model", "Body", "MODEL_VIEW")
+        rep = self.get_representation(body, rep_items)
+
+        bend_data = {"start_length": start_length, "end_length": end_length, "radius": radius, "angle": degrees(theta)}
+        return rep, bend_data