Uh60

.github/workflows/python-package-conda.yml

@@ -68,11 +68,11 @@ jobs:
           export MPLBACKEND=Agg
           python 2_sonata_IEA22.py
 
-      # - name: Run example 5
-      #   run: |
-      #     cd examples/5_NREL_5MW
-      #     export MPLBACKEND=Agg
-      #     python 5_sonata_NREL_5MW.py
+      - name: Run example 5
+        run: |
+          cd examples/5_UH60
+          export MPLBACKEND=Agg
+          python 5_uh60.py
 
       - name: Run example 6
         run: |

SONATA/cbm/classCBM.py

@@ -23,9 +23,14 @@
 from SONATA.cbm.cbm_utl import trsf_sixbysix
 from SONATA.cbm.classBeamSectionalProps import BeamSectionalProps
 from SONATA.cbm.classCBMConfig import CBMConfig
+
 from SONATA.cbm.display.display_mesh import plot_cells
 from SONATA.cbm.mesh.cell import Cell
 from SONATA.cbm.mesh.consolidate_mesh import consolidate_mesh_on_web
+from SONATA.cbm.mesh.mesh_core import gen_core_cells
+
+from SONATA.cbm.mesh.mesh_intersect import map_mesh_by_intersect_curve2d
+
 from SONATA.cbm.mesh.mesh_utils import (grab_nodes_of_cells_on_BSplineLst,
                                         sort_and_reassignID,)
 from SONATA.cbm.mesh.node import Node
@@ -343,7 +348,7 @@ def cbm_gen_mesh(self, **kwargs):
         global_minLen = round(self.refL / Resolution, 5)
 
         core_cell_area = 1.0 * global_minLen ** 2
-        # bw_cell_area = 0.7 * global_minLen ** 2
+        bw_cell_area = 0.7 * global_minLen ** 2
         web_consolidate_tol = 0.5 * global_minLen
 
         # ===================MESH SEGMENT
@@ -368,7 +373,7 @@ def cbm_gen_mesh(self, **kwargs):
             (web.wr_nodes, web.wr_cells) = grab_nodes_of_cells_on_BSplineLst(self.SegmentLst[web.ID+1].cells, web.BSplineLst)
 
             if not web.wl_nodes or not web.wl_cells or not web.wr_nodes or not web.wr_cells:  # in case there was no mesh in a segment
-                print('STATUS:\t No mesh on Web Interface ' + str(web.ID) + ' to be consolodated.')
+                print('STATUS:\t No mesh on Web Interface ' + str(web.ID) + ' to be consolidated.')
 
                 # This message gets printed in cases where the web doesn't get
                 # meshed because there is no isotropic layer.
@@ -386,6 +391,21 @@ def cbm_gen_mesh(self, **kwargs):
         for c in self.mesh:
             tmp.extend(c.split_quads())
         self.mesh = tmp
+        # ============= BALANCE WEIGHT - CUTTING HOLE ALGORITHM
+        if self.config.setup["BalanceWeight"]:
+            print("STATUS:\t Meshing Balance Weight")
+
+            self.mesh, boundary_nodes = map_mesh_by_intersect_curve2d(self.mesh, self.BW.Curve, self.BW.wire, global_minLen)
+            # boundary_nodes = merge_nodes_if_too_close(boundary_nodes,self.BW.Curve,global_minLen,tol=0.05))
+            [bw_cells, bw_nodes] = gen_core_cells(boundary_nodes, bw_cell_area)
+
+            for c in bw_cells:
+                c.structured = False
+                c.theta_3 = 0
+                c.MatID = self.config.bw["Material"]
+                c.calc_theta_1()
+
+            self.mesh.extend(bw_cells)
 
         # invert nodes list of all cell to make sure they are counterclockwise for vabs in the right coordinate system!
         for c in self.mesh:

SONATA/cbm/mesh/mesh_intersect.py

@@ -0,0 +1,245 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Mar 30 15:55:34 2017
+
+@author: TPflumm
+"""
+# Core Library modules
+import math
+
+# Third party modules
+# THE PYTHON SHAPELY MODULE:
+import shapely.geometry as shp_geom
+# PythonOCC Libraries
+from OCC.Core.BRepExtrema import BRepExtrema_DistShapeShape
+from OCC.Core.GCPnts import GCPnts_QuasiUniformAbscissa
+from OCC.Core.Geom2d import Geom2d_Line
+from OCC.Core.Geom2dAdaptor import Geom2dAdaptor_Curve
+from OCC.Core.Geom2dAPI import Geom2dAPI_InterCurveCurve
+from OCC.Core.gp import (gp_Dir2d, gp_Lin2d, gp_Pnt2d, gp_Vec2d,)
+
+# First party modules
+from SONATA.cbm.mesh.mesh_utils import \
+    remove_duplicates_from_list_preserving_order
+from SONATA.cbm.topo.BSplineLst_utils import findPnt_on_curve
+
+
+def map_node_on_curve(node, Curve2d, theta_11, distance=1e5, **kwargs):
+    """
+    maps a node onto a curve2d.
+
+    Parameters
+    ----------
+    node : TYPE
+        DESCRIPTION.
+    Curve2d : TYPE
+        DESCRIPTION.
+    theta_11 : TYPE
+        DESCRIPTION.
+    distance : TYPE, optional
+        DESCRIPTION. The default is 1e5.
+    **kwargs : TYPE
+        DESCRIPTION.
+
+    Returns
+    -------
+    None.
+
+    """
+
+    # ==================DIRECTION 1 ===============================
+    Line1 = Geom2d_Line(gp_Lin2d(node.Pnt2d, gp_Dir2d(math.cos(math.radians(theta_11)), math.sin(math.radians(theta_11)))))
+    # display.DisplayShape(Line1,color='BLACK')
+    Intersection1 = Geom2dAPI_InterCurveCurve(Curve2d, Line1)
+    dist = distance
+    for i in range(1, Intersection1.NbPoints() + 1):
+        if node.Pnt2d.Distance(Intersection1.Point(i)) < dist:
+            dist = node.Pnt2d.Distance(Intersection1.Point(i))
+            NewPnt1 = Intersection1.Point(i)
+
+    # ===================DIRECTION 2 ===============================
+    Line2 = Geom2d_Line(gp_Lin2d(node.Pnt2d, gp_Dir2d(-math.sin(math.radians(theta_11)), math.cos(math.radians(theta_11)))))
+    # display.DisplayShape(Line2,color='WHITE')
+    Intersection2 = Geom2dAPI_InterCurveCurve(Curve2d, Line2)
+    dist = distance
+    for i in range(1, Intersection2.NbPoints() + 1):
+        if node.Pnt2d.Distance(Intersection2.Point(i)) < dist:
+            dist = node.Pnt2d.Distance(Intersection2.Point(i))
+            NewPnt2 = Intersection2.Point(i)
+
+    # ==================CHOOSE NEW POINT===============================
+    try:
+        if node.Pnt2d.Distance(NewPnt1) > node.Pnt2d.Distance(NewPnt2):
+            node.Pnt2d = NewPnt2
+        else:
+            node.Pnt2d = NewPnt1
+
+    except NameError:
+        try:
+            node.Pnt2d = NewPnt2
+        except NameError:
+            node.Pnt2d = NewPnt1
+
+    return None
+
+
+def map_mesh_by_intersect_curve2d(mesh, curve2d, wire, global_minLen, **kwargs):
+    """
+    maps the mesh with and intersecting 2dcurve
+
+    Parameters
+    ----------
+    mesh : list of cells
+        DESCRIPTION.
+    curve2d : curve2d
+        DESCRIPTION.
+    wire : OCC.topods.wire
+        DESCRIPTION.
+    global_minLen : TYPE
+        DESCRIPTION.
+    **kwargs : TYPE
+        DESCRIPTION.
+
+    Returns
+    -------
+    None.
+
+    """
+    # KWARGS:
+    if kwargs.get("display") is not None:
+        display = kwargs.get("display")
+    else:
+        display = None
+
+    # # =========================INTERSECTING CELLS===============================
+    icells = []
+    for c in mesh:
+        Distance = BRepExtrema_DistShapeShape(wire, c.wire)
+        if Distance.Value() < 1e-5:
+            icells.append(c)
+
+    # Determine the interior nodes of Intersecting Cells
+    # ===================INTERIOR NODES of Intersecting Cells===================
+    Adaptor = Geom2dAdaptor_Curve(curve2d)
+    NbPoints = 60
+    discretization = GCPnts_QuasiUniformAbscissa(Adaptor, NbPoints)
+    shpPointLst2 = []
+    for j in range(1, NbPoints):
+        para = discretization.Parameter(j)
+        Pnt = gp_Pnt2d()
+        curve2d.D0(para, Pnt)
+        shpPointLst2.append((Pnt.X(), Pnt.Y()))
+
+    shp_Poly = shp_geom.Polygon(shpPointLst2)
+
+    inodes = []
+    for c in icells:
+        for n in c.nodes:
+            shpPnt = shp_geom.Point(n.coordinates[0], n.coordinates[1])
+            if shpPnt.within(shp_Poly):
+                inodes.append(n)
+                c.interior_nodes.append(n)
+
+    inodes = sorted(set(inodes), key=lambda Node: (Node.id))
+    #    for n in inodes:
+    #        display.DisplayShape(n.Pnt2d,color='Orange')
+
+    # =========================NONINTERSECTING INTERIOR CELLS===================
+    niicells = []
+    for c in mesh:
+        inodes = []
+        for n in c.nodes:
+            shpPnt = shp_geom.Point(n.coordinates[0], n.coordinates[1])
+            if shpPnt.within(shp_Poly):
+                inodes.append(n)
+        if len(inodes) == len(c.nodes):
+            niicells.append(c)
+
+    # for c in niicells:
+    # display.DisplayShape(c.wire,color='RED')
+
+    # =========================MAPPING==========================================
+    mapped_nodes = []
+    for c in icells:
+        if len(c.nodes) == 4 and len(c.interior_nodes) == 1:
+            node = c.interior_nodes[0]
+            if node not in mapped_nodes:
+                map_node_on_curve(node, curve2d, c.theta_11, display=display)
+                mapped_nodes.append(node)
+
+    for c in icells:
+        if len(c.nodes) == 3 and len(c.interior_nodes) == 1:
+            node = c.interior_nodes[0]
+            if node not in mapped_nodes:
+                map_node_on_curve(node, curve2d, c.theta_11, display=display)
+                mapped_nodes.append(node)
+
+    for c in icells:
+        if len(c.nodes) == 4 and len(c.interior_nodes) == 2:
+            node = c.interior_nodes[0]
+            if node not in mapped_nodes:
+                map_node_on_curve(node, curve2d, c.theta_11, display=display)
+                mapped_nodes.append(node)
+
+            node = c.interior_nodes[1]
+            if node not in mapped_nodes:
+                map_node_on_curve(node, curve2d, c.theta_11, display=display)
+                mapped_nodes.append(node)
+
+    rmcells = []
+    for c in icells:
+        if len(c.nodes) == 4 and len(c.interior_nodes) == 3:
+            s = set(c.nodes) - set(c.interior_nodes)
+            node = s.pop()
+            if node not in mapped_nodes:
+                map_node_on_curve(node, curve2d, c.theta_11, display=display)
+                mapped_nodes.append(node)
+            rmcells.append(c)
+
+    for c in icells:
+        if len(c.nodes) == 3 and len(c.interior_nodes) == 2:
+            t = set(c.nodes) - set(mapped_nodes)
+            if len(t) == 2:
+                s = set(c.nodes) - set(c.interior_nodes)
+                node = s.pop()
+
+                if node not in mapped_nodes:
+                    map_node_on_curve(node, curve2d, c.theta_11, display=display)
+                    mapped_nodes.append(node)
+                rmcells.append(c)
+
+    # =========================POPPING CELLS====================================
+    mesh = list(set(mesh) - set(niicells))
+    mesh = list(set(mesh) - set(rmcells))
+
+    rm2cells = []
+    for c in list(set(mesh) - set(rmcells)):
+        for n in c.nodes:
+            shpPnt = shp_geom.Point(n.coordinates[0], n.coordinates[1])
+            if shpPnt.within(shp_Poly):
+                rm2cells.append(c)
+
+    mesh = list(set(mesh) - set(rm2cells))
+
+    # TODO: If two nodes are mapped to the same point, then
+    # =========================SORTING MAPPED NODES=============================
+    uLst = []
+    mapped_nodes = list(set(mapped_nodes))
+    for n in mapped_nodes:
+        uLst.append(findPnt_on_curve(n.Pnt2d, curve2d))
+    mapped_nodes = [x for y, x in sorted(zip(uLst, mapped_nodes))]
+
+    # ====merge_mapped_nodes if to close:
+    tmp = []
+    for i, n1 in enumerate(mapped_nodes[0:], start=0):
+        n2 = mapped_nodes[i - 1]
+        v = gp_Vec2d(n1.Pnt2d, n2.Pnt2d)
+        magnitude = v.Magnitude()
+
+        if magnitude <= 0.001 * global_minLen:
+            n1 = n2
+        tmp.append(n1)
+
+    mapped_nodes = remove_duplicates_from_list_preserving_order(tmp)
+
+    return mesh, mapped_nodes

SONATA/cbm/topo/segment.py

@@ -4,9 +4,15 @@
 from SONATA.cbm.mesh.mesh_core import gen_core_cells
 from SONATA.cbm.mesh.mesh_utils import (grab_nodes_on_BSplineLst, remove_dublicate_nodes,
     remove_duplicates_from_list_preserving_order,)
+
+from OCC.Core.gp import gp_Vec2d
+
+
+
 from SONATA.cbm.topo.BSplineLst_utils import (copy_BSplineLst,
                                               get_BSplineLst_Pnt2d,
                                               reverse_BSplineLst,
+                                              get_BSplineLst_D2,
                                               trim_BSplineLst,)
 from SONATA.cbm.topo.layer import Layer
 from SONATA.cbm.topo.layer_utils import get_layer
@@ -170,6 +176,31 @@ def get_BsplineLst_plus(self, lid, SegmentLst, WebLst, layer_attr="a_BSplineLst"
 
         return (BSplineLst, start, end)
 
+    def det_weight_Pnt2d(self, s, t):
+        """
+        determine the position of the 2d Point that describes the balance
+        weight
+
+        Parameters
+        --------
+        s : float
+            nondimensional s position between Start S1 and End S2
+        t : float
+            distance in normaldirection left of the boundary_BesplineLst
+
+        Returns:
+        -------
+        p1 : gp_Pnt2d
+
+        """
+        p0, v1, v2 = get_BSplineLst_D2(self.BSplineLst, s, 0, 1)
+        # print('det_weight_Pnt2d:', p0.Coord(), 's', s)
+        v = gp_Vec2d(v1.Y(), -v1.X())
+        v.Normalize()
+        v.Multiply(t)
+        p1 = p0.Translated(v)
+        return p1
+
     def get_Pnt2d(self, S, SegmentLst, WebLst=None):
         """returns a Pnt2d for the coresponding layer number and the coordinate S"""
         a = self.final_Boundary_ivLst

examples/1_IEA15MW/1_sonata_IEA15.py

@@ -54,8 +54,8 @@
 
 # ===== User defined radial stations ===== #
 # Define the radial stations for cross sectional analysis (only used for flag_wt_ontology = True -> otherwise, sections from yaml file are used!)
-radial_stations =  [0., 0.01, 0.03, 0.05, 0.075, 0.15, 0.25, 0.3 , 0.4, 0.5 , 0.6 , 0.7 , 0.8 , 0.9 , 1.]
-# radial_stations = [.7]
+# radial_stations =  [0., 0.01, 0.03, 0.05, 0.075, 0.15, 0.25, 0.3 , 0.4, 0.5 , 0.6 , 0.7 , 0.8 , 0.9 , 1.]
+radial_stations = np.linspace(0., 1., 5, endpoint=True)
 # ===== Execute SONATA Blade Component Object ===== #
 # name          - job name of current task
 # filename      - string combining the defined folder directory and the job name

examples/2_IEA22MW/2_sonata_IEA22.py

@@ -58,7 +58,7 @@
 #               0.26530612, 0.28571429, 0.30612245, 0.32653061, 0.34693878, 0.36734694, 0.3877551 , 0.40816327, 0.42857143, 0.44897959, 0.46938776, 0.48979592, 0.51020408,
 #               0.53061224, 0.55102041, 0.57142857, 0.59183673, 0.6122449 , 0.63265306, 0.65306122, 0.67346939, 0.69387755, 0.71428571, 0.73469388, 0.75510204, 0.7755102 ,
 #               0.79591837, 0.81632653, 0.83673469, 0.85714286, 0.87755102, 0.89795918, 0.91836735, 0.93877551, 0.95918367, 1.]
-radial_stations = np.linspace(0., 1., 10, endpoint=True)
+radial_stations = np.linspace(0., 1., 5, endpoint=True)
 
 # ===== Execute SONATA Blade Component Object ===== #
 # name          - job name of current task