Merge pull request #12 from KIDEVLIN/SONATAmini

SONATAmini Updates

SONATA/cbm/classCBM.py

@@ -35,8 +35,13 @@
 from SONATA.cbm.topo.weight import Weight
 from SONATA.vabs.classStrain import Strain
 from SONATA.vabs.classStress import Stress
+from SONATA.cbm.topo.projection import (
+    chop_interval_from_layup, sort_layup_projection,)
 # from SONATA.vabs.classVABSConfig import VABSConfig
 
+from OCC.Core.Geom2dAPI import Geom2dAPI_InterCurveCurve
+from OCC.Core.gp import gp_Pnt2d
+
 try:
     import dolfin as do
     from SONATA.anbax.anbax_utl import build_dolfin_mesh, anbax_recovery, ComputeShearCenter, ComputeTensionCenter, ComputeMassCenter
@@ -178,6 +183,8 @@ def __init__(self, Configuration, materials=None, **kwargs):
         self.Theta = 0
         self.refL = get_BSplineLst_length(self.BoundaryBSplineLst)
 
+        self.cutoff_style = kwargs.get("cutoff_style")
+
     def _cbm_generate_SegmentLst(self, **kwargs):
         """
         psydo private method of the cbm class to generate the list of 
@@ -198,6 +205,56 @@ def _cbm_generate_SegmentLst(self, **kwargs):
         sorted(self.SegmentLst, key=getID)
         self.refL = get_BSplineLst_length(self.SegmentLst[0].BSplineLst)
         return None
+
+    def find_bspline_ends(self, bspline):
+        start_point = gp_Pnt2d()
+        end_point = gp_Pnt2d()
+        bspline.D0(bspline.FirstParameter(), start_point)
+        bspline.D0(bspline.LastParameter(), end_point)
+        return start_point, end_point
+
+    def check_bspline_intersections(self, Boundary_BSplineLst):
+        # Checking for bspline intersections in a bspline list not including start and end points
+        start_tol = 1e-5
+        intersection_points = []
+        intersected = False
+        for i in range(len(Boundary_BSplineLst)):
+            for j in range(len(Boundary_BSplineLst)):
+                if i != j:
+                    start1, end1 = self.find_bspline_ends(Boundary_BSplineLst[i])
+                    start2, end2 = self.find_bspline_ends(Boundary_BSplineLst[j])
+                    intersector = Geom2dAPI_InterCurveCurve(Boundary_BSplineLst[i],Boundary_BSplineLst[j])
+                    if intersector.NbPoints() > 0:
+                        for k in range(1, intersector.NbPoints()+1):
+                            if not intersector.Point(k).IsEqual(start1, start_tol) and not intersector.Point(k).IsEqual(end1, start_tol) and not intersector.Point(k).IsEqual(start2, start_tol) and not intersector.Point(k).IsEqual(end2, start_tol):
+                                intersected = True
+                                intersection_points.append(intersector.Point(k))
+        return [intersected, intersection_points]
+
+    def display_bsplinelst(self, bsplinelst, color = 'blue'):
+        for bspline in bsplinelst:
+            u_min, u_max = bspline.FirstParameter(), bspline.LastParameter()
+            # Extract points for plotting
+            num_points = 100  # Number of points to plot
+            u_values = [u_min + (u_max - u_min) * i / (num_points - 1) for i in range(num_points)]
+            x_values = [bspline.Value(u).X() for u in u_values]
+            y_values = [bspline.Value(u).Y() for u in u_values]
+            plt.plot(x_values, y_values, color = color)
+
+    def check_for_bspline_intersections(self, segment):
+        # Getting all boundary Bsplines
+        ivLst = chop_interval_from_layup(segment.boundary_ivLst, 0, 1)
+        ivLst = sort_layup_projection([ivLst])[0]
+        # Creating reference layer from the chopped ivLst
+        Boundary_BSplineLst = segment.ivLst_to_BSplineLst(ivLst)
+        [intersected, intersection_pnt] = self.check_bspline_intersections(Boundary_BSplineLst)
+        if intersected:
+            print("WARNING: There is an intersection in the structure.")
+            plt.figure()
+            self.display_bsplinelst(self.SegmentLst[0].BSplineLst, 'black')
+            self.display_bsplinelst(Boundary_BSplineLst, 'blue')
+            for points in intersection_pnt:
+                plt.plot(points.X(), points.Y(), 'x', color = 'red', linewidth = 4, markersize = 10)
 
     def cbm_gen_topo(self, **kwargs):
         """
@@ -212,8 +269,9 @@ def cbm_gen_topo(self, **kwargs):
         self._cbm_generate_SegmentLst(**kwargs)
         # Build Segment 0:
         self.SegmentLst[0].build_wire()
-        self.SegmentLst[0].build_layers(l0=self.refL, **kwargs)
+        self.SegmentLst[0].build_layers(l0=self.refL, cutoff_style = self.cutoff_style, **kwargs)
         self.SegmentLst[0].determine_final_boundary()
+        self.check_for_bspline_intersections(self.SegmentLst[0])
 
         # Build Webs:
         self.WebLst = []
@@ -229,8 +287,9 @@ def cbm_gen_topo(self, **kwargs):
                 seg.Segment0 = self.SegmentLst[0]
                 seg.WebLst = self.WebLst
                 seg.build_segment_boundary_from_WebLst(self.WebLst, self.SegmentLst[0])
-                seg.build_layers(self.WebLst, self.SegmentLst[0], l0=self.refL)
+                seg.build_layers(self.WebLst, self.SegmentLst[0], l0=self.refL, cutoff_style = self.cutoff_style)
                 seg.determine_final_boundary(self.WebLst, self.SegmentLst[0])
+                self.check_for_bspline_intersections(seg)
                 seg.build_wire()
 
         self.BW = None

SONATA/cbm/mesh/mesh_byprojection.py

@@ -5,11 +5,9 @@
 import numpy as np
 from OCC.Core.Geom2dAPI import (Geom2dAPI_PointsToBSpline,
                                 Geom2dAPI_ProjectPointOnCurve,)
-from OCC.Core.gp import gp_Pnt2d, gp_Vec2d, gp_Dir2d, gp_Lin2d
+from OCC.Core.gp import gp_Pnt2d, gp_Vec2d, gp_Dir2d
 from OCC.Core.Geom2d import Geom2d_Line
 from OCC.Core.Geom2dAPI import Geom2dAPI_InterCurveCurve
-from OCC.Core.Geom2dAPI import Geom2dAPI_Interpolate
-from OCC.Core.TColgp import TColgp_HArray1OfPnt2d
 
 # First party modules
 from SONATA.cbm.mesh.cell import Cell
@@ -35,7 +33,20 @@ def angle_between(v1, v2):
     return np.arccos(cos_theta)
 
 
-def projection_method(points, bspline, tolerance):
+def three_pnt_projection_method(points, bspline, min_dist):
+    """
+    Project node onto target Bspline along the bisecting vector based on the two adjacent nodes.
+    ----------
+    Input:
+    points: array of three points in order. The middle point is the point being projected and the
+    first and last points are the two points adjacent to the point being projected.
+    bspline: The targeted bspline.
+    tolerance: A distance tolerance. If the point is projected a distance farther than the tolerance,
+    it is not projected.
+
+    Output:
+    The point at which the bisecting vector intersects the target bspline.
+    """
     prev = points[0]
     curr = points[1]
     start_point_2d = gp_Pnt2d(curr[0],curr[1])
@@ -44,22 +55,23 @@ def projection_method(points, bspline, tolerance):
     vec2 = next - curr
     angle = angle_between(vec1,vec2)
 
+    # Find the two normal vectors pointing from the middle point to the adjacent points
     vec1_perp = np.array([-vec1[1], vec1[0]]) / np.linalg.norm([-vec1[1], vec1[0]])
     vec2_perp = np.array([-vec2[1], vec2[0]]) / np.linalg.norm([-vec2[1], vec2[0]])
 
+    # Finding the vector that bisects the two normal vectors. The project point will be in the direction of the bisecting vector.
     bisecting_vector = - np.sqrt((1 + np.dot(vec1_perp, vec2_perp)) / 2) * (vec1_perp + vec2_perp) / np.linalg.norm(vec1_perp + vec2_perp)
     direction_2d = gp_Dir2d(bisecting_vector[0],bisecting_vector[1])
     line_2d = Geom2d_Line(start_point_2d, direction_2d)
+
+    # Finding where the bisecting vector intersects the target Bspline
     intersector = Geom2dAPI_InterCurveCurve(line_2d, bspline)
-    # line_2d_neg = Geom2d_Line(start_point_2d, direction_2d)
-    # if not intersector:
-    #     intersector = Geom2dAPI_InterCurveCurve(line_2d_neg, bspline)
     if intersector.NbPoints() > 0:
         intersection_point = intersector.Point(1)
         distance = np.sqrt((curr[0]-intersection_point.X())**2+(curr[1]-intersection_point.Y())**2)
-        if distance < tolerance*(2+np.cos(angle)):
-            return intersection_point
-    return None
+        if distance < min_dist:
+            return [intersection_point, distance, True]
+    return [None, None, False]
 
 def mesh_by_projecting_nodes_on_BSplineLst(a_BSplineLst, a_nodes, b_BSplineLst, layer_thickness, tol=1e-2, crit_angle=95, LayerID=0, refL=1.0, **kw):
     """
@@ -119,6 +131,9 @@ def mesh_by_projecting_nodes_on_BSplineLst(a_BSplineLst, a_nodes, b_BSplineLst,
         pPara = []
         pIdx = []
         projected = False
+
+        # Projects current node onto the target Bspline where the vector between the current node
+        # and the projected node is perpendicular.
         for idx, item in enumerate(b_BSplineLst):
             first = item.FirstParameter()
             last = item.LastParameter()
@@ -135,21 +150,27 @@ def mesh_by_projecting_nodes_on_BSplineLst(a_BSplineLst, a_nodes, b_BSplineLst,
                     projected = True
                 else:
                     None
+        # If no node is found that creates a vector perpendicular to the target Bspline, a different
+        # projection method is tried using the current node and the neighboring nodes to create a
+        # bisecting vector. The point at which the bisecting vector intersects the target Bspline
+        # is where the projected point is placed. This helps for projecting nodes on corners.
         if not projected and i-1 != 0 and not i-1 > len(prj_nodes)-2:
+            min_distance = np.inf
             for idx, item in enumerate(b_BSplineLst):
                 prev_point = np.array([prj_nodes[i-2].coordinates[0], prj_nodes[i-2].coordinates[1]])
                 curr_point = np.array([prj_nodes[i-1].coordinates[0], prj_nodes[i-1].coordinates[1]])
                 next_point = np.array([prj_nodes[i].coordinates[0], prj_nodes[i].coordinates[1]])
-                projected_point = projection_method([prev_point,curr_point,next_point], item, distance*10)
-                if projected_point:
-                    pPnts.append(projected_point)
-                    pPara.append(node.parameters[2])
-                    pIdx.append(idx)
-                    plt.plot(projected_point.X(),projected_point.Y(), color = 'purple', marker = 'x', markersize = 12)
-                    node.corner = False
-                    break
-        for point in pPnts:
-            plt.plot(point.X(),point.Y(), color = 'blue', marker = 'o', markersize = 10)
+                [tmp_projected_point, tmp_min_distance, projected] = three_pnt_projection_method([prev_point,curr_point,next_point], item, min_distance)
+                if projected:
+                    projected_point = tmp_projected_point
+                    min_distance = tmp_min_distance
+                    tmp_idx = idx
+
+            if projected_point:
+                pPnts.append(projected_point)
+                pPara.append(node.parameters[2])
+                pIdx.append(tmp_idx)
+                node.corner = False
 
 
         # ==================making sure the pPnts are unique:
@@ -554,9 +575,5 @@ def mesh_by_projecting_nodes_on_BSplineLst(a_BSplineLst, a_nodes, b_BSplineLst,
 
         display.View_Top()
         display.FitAll()
-    # for node in a_nodes:
-    #          plt.plot(node.coordinates[0],node.coordinates[1], color = 'blue', marker = 'o', markersize = 10)
-    # for node in b_nodes:
-    #         plt.plot(node.coordinates[0],node.coordinates[1], color = 'green', marker = 'o', markersize = 10)
 
     return a_nodes, b_nodes, cellLst

SONATA/cbm/topo/layer.py

@@ -109,72 +109,12 @@ def __setstate__(self, state):
     def build_wire(self):  # Builds TopoDS_Wire from connecting BSplineSegments and returns it
         self.wire = build_wire_from_BSplineLst(self.BSplineLst)
 
-    def display_bsplinelst(self, bsplinelst, color = 'red'):
-        for bspline in bsplinelst:
-            u_min, u_max = bspline.FirstParameter(), bspline.LastParameter()
-            # Extract points for plotting
-            num_points = 100  # Number of points to plot
-            u_values = [u_min + (u_max - u_min) * i / (num_points - 1) for i in range(num_points)]
-            x_values = [bspline.Value(u).X() for u in u_values]
-            y_values = [bspline.Value(u).Y() for u in u_values]
-
-            plt.plot(x_values, y_values, color = color)
-
-    # Function to check if two line segments (p1, q1) and (p2, q2) intersect
-    def do_intersect(self, p1, q1, p2, q2):
-        def orientation(p, q, r):
-            val = (float(q[1] - p[1]) * (r[0] - q[0])) - (float(q[0] - p[0]) * (r[1] - q[1]))
-            if val > 0:
-                return 1  # Clockwise
-            elif val < 0:
-                return 2  # Counterclockwise
-            else:
-                return 0  # Collinear
-
-        def on_segment(p, q, r):
-            if min(p[0], q[0]) <= r[0] <= max(p[0], q[0]) and min(p[1], q[1]) <= r[1] <= max(p[1], q[1]):
-                return True
-            return False
-
-        o1 = orientation(p1, q1, p2)
-        o2 = orientation(p1, q1, q2)
-        o3 = orientation(p2, q2, p1)
-        o4 = orientation(p2, q2, q1)
-
-        # General case
-        if o1 != o2 and o3 != o4:
-            return True
-
-        # Special cases
-        if o1 == 0 and on_segment(p1, q1, p2):
-            return True
-        if o2 == 0 and on_segment(p1, q1, q2):
-            return True
-        if o3 == 0 and on_segment(p2, q2, p1):
-            return True
-        if o4 == 0 and on_segment(p2, q2, q1):
-            return True
-
-        return False
-
-    # Function to check if a shape intersects itself
-    def shape_intersects_itself(self, coords):
-        num_points = len(coords)
-        for i in range(num_points - 1):  # Loop to (num_points - 1) to avoid out-of-bounds access
-            for j in range(i + 2, num_points - (1 if i == 0 else 0)):  # Adjusting the range to avoid out-of-bounds
-                if self.do_intersect(coords[i], coords[i + 1], coords[j], coords[(j + 1) % num_points]):
-                    return True
-        return False 
-
     def build_layer(self, l0=1):
         npArray = discretize_BSplineLst(self.Boundary_BSplineLst, 1.2e-6 * l0)
         self.offlinepts = shp_parallel_offset(npArray, self.thickness, self.join_style)
-        if self.shape_intersects_itself(self.offlinepts[1:-2]):
-            print("WARNING: THERE IS AN INTERSECTION IN THE STRUCTURE")
         OffsetBSplineLst = BSplineLst_from_dct(self.offlinepts, angular_deflection=15, tol_interp=1e-8 * l0, cutoff_style = 2)
         OffsetBSplineLst = cutoff_layer(self.Boundary_BSplineLst, OffsetBSplineLst, self.S1, self.S2, self.cutoff_style)
         self.BSplineLst = OffsetBSplineLst
-        self.display_bsplinelst(self.BSplineLst)
 
     def determine_a_nodes(self, SegmentLst, global_minLen, display=None):
         """ """

SONATA/cbm/topo/offset.py

@@ -8,11 +8,70 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import shapely.geometry as shp
+from shapely.geometry import MultiPolygon
 
 # Local modules
 from .utils import (P2Pdistance, Polygon_orientation,
                     calc_DCT_angles, isclose, unique_rows,)
 
+def get_largest_polygon(multipolygon):
+    # Ensure the input is a MultiPolygon
+    if not isinstance(multipolygon, MultiPolygon):
+        raise ValueError("Input must be a MultiPolygon object")
+
+    # Initialize variables to track the largest polygon
+    largest_polygon = None
+    largest_area = 0
+
+    # Iterate through each polygon in the MultiPolygon
+    for polygon in multipolygon.geoms:
+        # Compute the area of the current polygon
+        area = polygon.area
+        # Check if this polygon has a larger area than the current largest
+        if area > largest_area:
+            largest_area = area
+            largest_polygon = polygon
+
+    return largest_polygon
+
+def combine_close_points(points, tolerance, length_threshold):
+    def distance(p1, p2):
+        return np.linalg.norm(p1 - p2)
+
+    def segment_length(start_idx, end_idx):
+        # Calculate the length of the segment from points[start_idx] to points[end_idx]
+        segment_length = 0
+        for i in range(start_idx, end_idx):
+            segment_length += distance(points[i], points[i+1])
+        return segment_length
+
+    combined_points = []
+    i = 0
+
+    while i < len(points):
+        close_points = [points[i]]
+        j = i + 1
+
+        while j < len(points) and distance(points[j], points[j-1]) <= tolerance:
+            close_points.append(points[j])
+            j += 1
+
+        if len(close_points) > 1:
+            # Check if the length of the segment is within the threshold
+            segment_len = segment_length(i, j-1)
+            if segment_len < length_threshold and segment_len > 3 * tolerance:
+                # Keep the middle point if the length is less than the threshold
+                middle_index = len(close_points) // 2
+                combined_points.append(close_points[middle_index])
+            else:
+                # Keep all points if the length is greater than or equal to the threshold
+                combined_points.extend(close_points)
+        else:
+            combined_points.append(points[i])
+
+        i = j
+
+    return np.array(combined_points)
 
 # Function to check if two line segments (p1, q1) and (p2, q2) intersect
 def do_intersect(p1, q1, p2, q2):
@@ -72,12 +131,10 @@ def shp_parallel_offset(arrPts, dist, join_style=1, side="right", res=16):
             afpoly = shp.Polygon(arrPts)
             noffafpoly = afpoly.buffer(-dist)  # Inward offset
             data = np.array(noffafpoly.exterior.xy).T
-        except:
-            intersects = True
-            while intersects:
-                [intersects, arrPts] = shape_intersects_itself(arrPts)
+        except AttributeError as e:
             afpoly = shp.Polygon(arrPts)
-            noffafpoly = afpoly.buffer(-dist)  # Inward offset
+            noffaf_multipoly = afpoly.buffer(-dist)  # Inward offset
+            noffafpoly = get_largest_polygon(noffaf_multipoly) # If there are overlaps find polygon with largest area
             data = np.array(noffafpoly.exterior.xy).T
 
     else:
@@ -152,9 +209,8 @@ def shp_parallel_offset(arrPts, dist, join_style=1, side="right", res=16):
             Refinement = True
         else:
             Refinement = False
-        # if closed:
-        #     np.vstack((data, data[0]))
-
+    # Combine close points around corners
+    data = combine_close_points(data, dist**2*100, dist*4)
     return data