Select target elements using Point in Polygon Ray Casting Algorithm.

Ayush Kumar
Ayush Kumar Member, Moderator, Employee Posts: 450
100 Answers 250 Likes 100 Comments First Anniversary
✭✭✭✭
edited September 2023 in Structures

Components 1, 2 and 3 are connected to the base plate as bonded contacts. How can I select the element faces on the base plate which are attached to the component base?

Geometry

Output

Comments

  • Ayush Kumar
    Ayush Kumar Member, Moderator, Employee Posts: 450
    100 Answers 250 Likes 100 Comments First Anniversary
    ✭✭✭✭
    edited September 2023

    The following script uses Point in Polygon Ray Casting Algorithm to search for nodes of the component boundary in the element faces on the base plate [(Method-I)] and vice-versa [(Method-II)]. The output is a named selection which contains common elements of [Method-I] and [Method-II].

    Prerequisites:

    1. The components must be connected to the base plate via a bonded connection, with component faces as "Target"
    2. The base plate face must be scoped in a named selection called "base_face"
    3. Script considers a projection in X-Y plane.

    Please Note: This script must be used only as a reference script and users must adapt it to their own requirements. This script has not undergone robust testing, so there could be bugs. Please use at your own risk.

    import math
    from collections import defaultdict
    
    INTERIM_NS = []
    ALL_ELEMENT_FACES = []
    
    contact_face_type = "Target"
    base_face = ExtAPI.DataModel.GetObjectsByName('base_face')[0]
    mesh = ExtAPI.DataModel.MeshDataByName("Global")
    all_base_nodes_ns = base_face.CreateNodalNamedSelection()
    INTERIM_NS.append(all_base_nodes_ns)
    all_base_nodes_ns.Name = "all_base_nodes"
    all_base_nodes = all_base_nodes_ns.Ids
    
    def convert_nodes_to_element_faces(nodal_ns):
        ns_efs = Model.AddNamedSelection()
        ns_efs.ScopingMethod = GeometryDefineByType.Worksheet
    
        ns_efs.GenerationCriteria.Add(None)
        ns_efs.GenerationCriteria[0].EntityType = SelectionType.MeshNode
        ns_efs.GenerationCriteria[0].Criterion = SelectionCriterionType.NamedSelection
        ns_efs.GenerationCriteria[0].Operator = SelectionOperatorType.Equal
        ns_efs.GenerationCriteria[0].Value = nodal_ns
    
        ns_efs.GenerationCriteria.Add(None)
        ns_efs.GenerationCriteria[1].Action = SelectionActionType.Convert
        ns_efs.GenerationCriteria[1].EntityType = SelectionType.MeshElementFace
        ns_efs.Generate()
        return ns_efs
    
    def convert_faces_to_mesh_nodes(faces_ns):
        ns_component = Model.AddNamedSelection()
        ns_component.ScopingMethod = GeometryDefineByType.Worksheet
        ns_component.GenerationCriteria.Add(None)
        ns_component.GenerationCriteria[0].EntityType = SelectionType.GeoFace
        ns_component.GenerationCriteria[0].Criterion = SelectionCriterionType.NamedSelection
        ns_component.GenerationCriteria[0].Operator = SelectionOperatorType.Equal
        ns_component.GenerationCriteria[0].Value = faces_ns
    
        ns_component.GenerationCriteria.Add(None)
        ns_component.GenerationCriteria[1].Action = SelectionActionType.Convert
        ns_component.GenerationCriteria[1].EntityType = SelectionType.GeoEdge
    
        ns_component.GenerationCriteria.Add(None)
        ns_component.GenerationCriteria[2].Action = SelectionActionType.Convert
        ns_component.GenerationCriteria[2].EntityType = SelectionType.MeshNode
        ns_component.Generate()
        return ns_component
    
    def convert_to_element_faces(nodal_ns, from_entity=SelectionType.MeshNode):
        ns_efs = Model.AddNamedSelection()
        ns_efs.ScopingMethod = GeometryDefineByType.Worksheet
    
        ns_efs.GenerationCriteria.Add(None)
        ns_efs.GenerationCriteria[0].EntityType = from_entity
        ns_efs.GenerationCriteria[0].Criterion = SelectionCriterionType.NamedSelection
        ns_efs.GenerationCriteria[0].Operator = SelectionOperatorType.Equal
        ns_efs.GenerationCriteria[0].Value = nodal_ns
    
        ns_efs.GenerationCriteria.Add(None)
        ns_efs.GenerationCriteria[1].Action = SelectionActionType.Convert
        ns_efs.GenerationCriteria[1].EntityType = SelectionType.MeshElementFace
        ns_efs.Generate()
        return ns_efs
    
    def convert_selection_to_polygon(selection):
        element_ids = selection.Ids
        element_face_indices = selection.ElementFaceIndices
        ElementFaceWrapper = Ansys.ACT.Common.Mesh.ElementFaceWrapper
        polygon = defaultdict(list)
        for id, index in zip(element_ids, element_face_indices):
            el_face = ElementFaceWrapper(mesh, mesh.ElementById(id), index)
            polygon[id] = ([Point(mesh.NodeById(nid).X, mesh.NodeById(nid).Y) for nid in el_face.NodeIds], index)
        return polygon
    
    """ Sort random coordinates to form a closed polygon """
    
    def find_centroid(points):
        # Calculate the centroid of a list of (x, y) points
        num_points = len(points)
        sum_x = sum(point.x for point in points)
        sum_y = sum(point.y for point in points)
        centroid_x = sum_x / num_points
        centroid_y = sum_y / num_points
        return (centroid_x, centroid_y)
    
    def calculate_polar_angle(point, centroid):
        # Calculate the polar angle of a point with respect to the centroid
        delta_x = point.x - centroid[0]
        delta_y = point.y - centroid[1]
        return math.atan2(delta_y, delta_x)
    
    def sort_points_to_form_polygon(points):
        centroid = find_centroid(points)
        points.sort(key=lambda point: calculate_polar_angle(point, centroid))
        return points
    
    class Point:
        def __init__(self, x, y):
            self.x = x
            self.y = y
    
    class Line:
        def __init__(self, p1, p2):
            self.p1 = p1
            self.p2 = p2
    
    def cross_product(p1, p2, p3):
        return (p2.x - p1.x) * (p3.y - p1.y) - (p2.y - p1.y) * (p3.x - p1.x)
    
    def is_on_line(side, point):
        m = (side.p1.y - side.p2.y) / (side.p1.x - side.p2.x)
        if (point.y - side.p1.y) == m * (point.x - side.p1.x):
            return True
        return False
    
    def direction(p1, p2, p3):
        cp = cross_product(p1, p2, p3)
        if cp == 0:
            # Collinear
            return 0
        elif cp < 0:
            # Anti-clockwise / Left side
            return 2
        # Clockwise / Right side
        return 1
    
    def intersects(l1, l2):
        dir1 = direction(l1.p1, l1.p2, l2.p1)
        dir2 = direction(l1.p1, l1.p2, l2.p2)
        dir3 = direction(l2.p1, l2.p2, l1.p1)
        dir4 = direction(l2.p1, l2.p2, l1.p2)
    
        # When intersecting
        if dir1 != dir2 and dir3 != dir4:
            return True
        # When p1 of line2 is on line1
        if dir1 == 0 and is_on_line(l1, l2.p1):
            return True
        # When p2 of line2 is on line1
        if dir2 == 0 and is_on_line(l1, l2.p2):
            return True
        # When p1 of line1 is on line2
        if dir3 == 0 and is_on_line(l2, l1.p1):
            return True
        # When p2 of line1 is on line2
        if dir4 == 0 and is_on_line(l2, l1.p2):
            return True
        return False
    
    """ Method to check if a point lies within a polygon using ``Ray Casting Algorithm`` """
    
    def is_inside_polygon(point, polygon):
        segments = len(polygon)
        ray = Line(point, Point(1.0e12, point.y))
        n_intersects = 0
        for index in range(segments):
            v1 = index
            v2 = 0 if (index == (segments - 1)) else (index + 1)
            side = Line(polygon[v1], polygon[v2])
            if intersects(side, ray):
                if direction(side.p1, side.p2, point) == 0:
                    return is_on_line(side, point)
                n_intersects += 1
    
        if n_intersects % 2 == 0:
            return False
        return True
    
    
  • Ayush Kumar
    Ayush Kumar Member, Moderator, Employee Posts: 450
    100 Answers 250 Likes 100 Comments First Anniversary
    ✭✭✭✭
    edited September 2023

    Please Note: This script must be used only as a reference script and users must adapt it to their own requirements. This script has not undergone robust testing, so there could be bugs. Please use at your own risk.

    Part - II:

    def main():
        contacts = ExtAPI.DataModel.GetObjectsByType(Ansys.ACT.Automation.Mechanical.Connections.ContactRegion)
        bonded_contacts = filter(lambda x: (x.ContactType == ContactType.Bonded), contacts)
    
        for index, bonded_contact in enumerate(bonded_contacts, start=1):
            if bonded_contact.Suppressed == True:
                continue
            ns_component_face = Model.AddNamedSelection()
            ns_component_face.Location = bonded_contact.SourceLocation if contact_face_type == "Source" else bonded_contact.TargetLocation
            ns_component_face.Name = "face_component_%s" % index
            INTERIM_NS.append(ns_component_face)
    
            """ Create Polygon Loop for base nodes in component loop """
            ns_component = convert_faces_to_mesh_nodes(ns_component_face)
            ns_component.Name = "loop_component_%s" % index
            INTERIM_NS.append(ns_component)
    
            """ Rough filter all base nodes to nodes near the component """
            x_min = min([mesh.NodeById(nid).X for nid in ns_component.Ids])
            x_max = max([mesh.NodeById(nid).X for nid in ns_component.Ids])
            y_min = min([mesh.NodeById(nid).Y for nid in ns_component.Ids])
            y_max = max([mesh.NodeById(nid).Y for nid in ns_component.Ids])
            temp_filt_nodes = filter(
                lambda x: (mesh.NodeById(x).X >= x_min and mesh.NodeById(x).X <= x_max and mesh.NodeById(
                    x).Y >= y_min and mesh.NodeById(x).Y <= y_max), all_base_nodes)
    
            """ Create and sort polygon with coordinates """
            random_polygon_points = [Point(mesh.NodeById(nid).X, mesh.NodeById(nid).Y) for nid in ns_component.Ids]
            polygon = sort_points_to_form_polygon(random_polygon_points)
    
            """ Method I - check for base nodes in component polygon """
            nodes_in_polygon = []
            for nid in temp_filt_nodes:
                point = Point(mesh.NodeById(nid).X, mesh.NodeById(nid).Y)
                if is_inside_polygon(point, polygon):
                    nodes_in_polygon.append(nid)
    
            sm = ExtAPI.SelectionManager.CreateSelectionInfo(SelectionTypeEnum.MeshNodes)
            sm.Ids = nodes_in_polygon
    
            ns_nodes = Model.AddNamedSelection()
            ns_nodes.Location = sm
            ns_nodes.Name = "base_face_nodes_component_%s" % index
            INTERIM_NS.append(ns_nodes)
    
            ns_efs_filt_1 = convert_nodes_to_element_faces(ns_nodes)
            ns_efs_filt_1.Name = "efs_method_1_component_%s" % index
            INTERIM_NS.append(ns_efs_filt_1)
    
            """ Method II - checking for component boundary nodes in base polygon """
    
            temp_sm = ExtAPI.SelectionManager.CreateSelectionInfo(SelectionTypeEnum.MeshNodes)
            temp_sm.Ids = temp_filt_nodes
            temp_ns = Model.AddNamedSelection()
            temp_ns.Location = temp_sm
            temp_ns.Name = "filtered_nodes_component_%s" % index
            INTERIM_NS.append(temp_ns)
    
            ns_temp = convert_nodes_to_element_faces(temp_ns)
            ns_temp.Name = "filtered_element_faces_component_%s" % index
            INTERIM_NS.append(ns_temp)
    
            base_polygons = convert_selection_to_polygon(ns_temp)
    
            filter_polygon = defaultdict(list)
            for nid in ns_component.Ids:
                point = Point(mesh.NodeById(nid).X, mesh.NodeById(nid).Y)
                p_count = 0
                while p_count < len(base_polygons):
                    if is_inside_polygon(point, base_polygons.values()[p_count][0]):
                        eid = base_polygons.keys()[p_count]
                        filter_polygon[eid] = base_polygons[eid]
                        base_polygons.pop(eid)
                        break
                    p_count += 1
    
            sm = ExtAPI.SelectionManager.CreateSelectionInfo(SelectionTypeEnum.MeshElements)
            sm.Ids = filter_polygon.keys()
    
            ns_elem = Model.AddNamedSelection()
            ns_elem.Location = sm
            ns_elem.Name = "filt_elements_component_%s" % index
            INTERIM_NS.append(ns_elem)
    
            ef_indices = [ef_index[1] for ef_index in filter_polygon.values()]
    
            sm_ef_indices = ExtAPI.SelectionManager.CreateSelectionInfo(SelectionTypeEnum.MeshElementFaces)
            sm_ef_indices.Ids = filter_polygon.keys()
            sm_ef_indices.ElementFaceIndices = ef_indices
            ns_efs_filt_2 = Model.AddNamedSelection()
            ns_efs_filt_2.Location = sm_ef_indices
            ns_efs_filt_2.Name = "efs_method_2_component__%s" % index
            INTERIM_NS.append(ns_efs_filt_2)
    
            """ Combine selected element faces from both the Methods """
            combined_ns = Model.AddNamedSelection()
            combined_ns.ScopingMethod = GeometryDefineByType.Worksheet
            combined_ns.GenerationCriteria.Add(None)
            combined_ns.GenerationCriteria[0].EntityType = SelectionType.MeshElementFace
            combined_ns.GenerationCriteria[0].Criterion = SelectionCriterionType.NamedSelection
            combined_ns.GenerationCriteria[0].Operator = SelectionOperatorType.Equal
            combined_ns.GenerationCriteria[0].Value = ns_efs_filt_1
    
            combined_ns.GenerationCriteria.Add(None)
            combined_ns.GenerationCriteria[1].EntityType = SelectionType.MeshElementFace
            combined_ns.GenerationCriteria[1].Criterion = SelectionCriterionType.NamedSelection
            combined_ns.GenerationCriteria[1].Operator = SelectionOperatorType.Equal
            combined_ns.GenerationCriteria[1].Value = ns_efs_filt_2
            combined_ns.Generate()
            combined_ns.Name = "element_faces_component_%s" % index
            ALL_ELEMENT_FACES.append(combined_ns)
    
        group = Tree.Group(INTERIM_NS)
        group.Name = "Temp_NS"
    
        group = Tree.Group(ALL_ELEMENT_FACES)
        group.Name = "ELEMENT_FACES"
    
        Model.NamedSelections.Activate()
    
    
    main()