Stress Example, Bug Fixes, Additional Options

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

@@ -65,3 +65,8 @@ jobs:
         run: |
           cd examples/5_IEA_5MW
           python 5_sonata_IEA5.py
+
+      - name: Run example 6
+        run: |
+          cd examples/6_beam_stress
+          python 6_sonata_beam_stress.py

README.md

@@ -13,13 +13,27 @@ SONATA can be run on mac and linux machines. No Windows installation is supporte
 
 At NREL (and possibly at other institutes), first disconnect from vpn client during installation in order to avoid remote server error when trying to retrieve URLs for installation.
 
-First setup an anaconda environment, here named sonata-env, activate it, and add the pythonocc library (v7.4.1)
+First setup an anaconda environment, here named sonata-env, activate it, and add the pythonocc library (v7.4.1).
+You should do this in the folder that you want to clone SONATA to.
 
 ```
 conda config --add channels conda-forge
 conda config --add channels tpaviot
-conda env create --name sonata-env -f https://raw.githubusercontent.com/ptrbortolotti/SONATA/master/environment.yaml python=3.9
+git clone [email protected]:ptrbortolotti/SONATA.git
+cd SONATA
+```
+If you have a mac with a newer chip, run the following:
+```
+CONDA_SUBDIR=osx-64 conda env create --name sonata-env -f environment.yaml
+conda activate sonata-env
+conda config --env --set subdir osx-64 # run with sonata-env active.
+cd ..
+```
+Otherwise, you should be able to run:
+```
+conda env create --name sonata-env -f environment.yaml
 conda activate sonata-env
+cd ..
 ```
 
 Next, download the solvers VABS (commercial, use wine to run it on mac/linux systems) or in the same conda environment compile ANBA4 (open-source)
@@ -31,10 +45,9 @@ pip install -e .
 cd ..
 ```
 
-Finally, go to the folder where you want to clone SONATA and type:
+Go back to where you cloned SONATA and type:
 
 ```
-git clone [email protected]:ptrbortolotti/SONATA.git
 cd SONATA
 pip install -e .
 ```
@@ -63,4 +76,4 @@ python 1_sonata_IEA15.py
 
 **Pflumm, T., Garre, W., Hajek, M.:** A Preprocessor for Parametric Composite Rotor Blade Cross-Sections, 44th European Rotorcraft Forum, Delft, The Netherlands, 2018  [[pdf]](docs/Pflumm,%20T.%20-%20A%20Preprocessor%20for%20Parametric%20Composite%20Rotor%20Blade%20Cross-Sections%20(2018,%20ERF).pdf) [[more…\]](https://mediatum.ub.tum.de/604993?query=Pflumm&show_id=1455385) [[BibTeX\]](https://mediatum.ub.tum.de/export/1455385/bibtex)
 
-**Pflumm, T., Rex, W., Hajek, M.:** Propagation of Material and Manufacturing Uncertainties in Composite Helicopter Rotor Blades, 45th European Rotorcraft Forum, Warsaw, Poland, 2019 [[more…\]](https://mediatum.ub.tum.de/1520025) [BibTeX\]](https://mediatum.ub.tum.de/export/1520025/bibtex)
+**Pflumm, T., Rex, W., Hajek, M.:** Propagation of Material and Manufacturing Uncertainties in Composite Helicopter Rotor Blades, 45th European Rotorcraft Forum, Warsaw, Poland, 2019 [[more…\]](https://mediatum.ub.tum.de/1520025) [BibTeX\]](https://mediatum.ub.tum.de/export/1520025/bibtex)

SONATA/cbm/classCBM.py

@@ -413,7 +413,7 @@ def cbm_run_anbax(self):
         except:
             print('\n')
             print('==========================================\n\n')
-            print('Error, Anba4 wrapper called, but ')
+            print('Error, Anba4 wrapper called, likely ')
             print('Anba4 _or_ Dolfin are not installed\n\n')
             print('==========================================\n\n')
 

SONATA/cbm/display/display_mesh.py

@@ -55,11 +55,12 @@ def plot_mesh(nodes, elements, theta_11, data, data_name, materials, title=None,
         cmap_name, colors, N=6)
 
     elif data_name == 'MatID':
-        cmap = a=plt.cm.get_cmap()
+        cmap = plt.cm.get_cmap()
         # extract all colors from the .jet map
         cmaplist = [cmap(i) for i in range(cmap.N)]
         cmap = LinearSegmentedColormap.from_list('Custom cmap', cmaplist, max(data))
 
+        cmap = cm.get_cmap('tab20',max(data))
     else:
         cmap = plt.cm.get_cmap()
 
@@ -72,9 +73,6 @@ def plot_mesh(nodes, elements, theta_11, data, data_name, materials, title=None,
         vmax = kw["vmax"]
     else:
         vmax = None
-    from palettable.cartocolors.qualitative import Prism_9
-    from matplotlib.colors import ListedColormap
-    cmap = cm.get_cmap('tab20',max(data))
 
     fig, ax = plt.subplots(1,1,figsize=(7.5, 6))
     patches = []
@@ -95,7 +93,7 @@ def plot_mesh(nodes, elements, theta_11, data, data_name, materials, title=None,
     p.set_clim(vmin, vmax)
     _ = ax.add_collection(p)
 
-    cbar = fig.colorbar(p, ax=ax, drawedges=True)
+    cbar = fig.colorbar(p, ax=ax, drawedges=(data_name == 'MatID'))
     cbar.ax.set_ylabel(data_name)
 
     if data_name == "MatID":

SONATA/cbm/mesh/mesh_byprojection.py

@@ -13,7 +13,8 @@
 from SONATA.cbm.mesh.cell import Cell
 from SONATA.cbm.mesh.node import Node
 from SONATA.cbm.topo.BSplineLst_utils import (
-    ProjectPointOnBSplineLst, find_BSplineLst_coordinate, intersect_BSplineLst_with_BSpline,)
+    ProjectPointOnBSplineLst, find_BSplineLst_coordinate,
+    intersect_BSplineLst_with_BSpline,discretize_BSplineLst)
 from SONATA.cbm.topo.utils import point2d_list_to_TColgp_Array1OfPnt2d
 
 def display_bsplinelst(bsplinelst, color = 'red'):
@@ -73,7 +74,10 @@ def three_pnt_projection_method(points, bspline, 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):
+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):
     """
     *function to mesh the SONATA topologies by projecting nodes onto the 
     generated BSplineLists.   
@@ -150,16 +154,33 @@ 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:
+        if not projected and (closed_a or (i-1 != 0 and i-1 <= len(prj_nodes)-2)):
+            # Requirements: not projected by previous method
+            # Section must be closed or need to have a point on both sides
             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]])
+
+                # If section is not closed, then only get here if i-2 >= 0
+                # if section is closed i-2 can be -1, which is the correct
+                # wrap around.
+                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]])
+
+                if closed_a and i-1 > len(prj_nodes)-2:
+                    # Since section is closed wrap around the index
+                    next_point = np.array([prj_nodes[0].coordinates[0],
+                                           prj_nodes[0].coordinates[1]])
+                else:
+                    next_point = np.array([prj_nodes[i].coordinates[0],
+                                           prj_nodes[i].coordinates[1]])
+
                 [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
@@ -172,6 +193,42 @@ def mesh_by_projecting_nodes_on_BSplineLst(a_BSplineLst, a_nodes, b_BSplineLst,
                 pIdx.append(tmp_idx)
                 node.corner = False
 
+        elif not projected and not closed_a:
+            print("WARNING: Doing discrete nearest point for projection."
+                  + " First or last point failed to project.")
+            # It is better to do this than just skip the projected node
+            # from looking at a few cases.
+
+            discrete_pts = len(b_BSplineLst) * [None]
+            discrete_pts_dist = len(b_BSplineLst) * [None]
+
+            for idx, item in enumerate(b_BSplineLst):
+                discrete_curve = discretize_BSplineLst([item], 1.0e-7*refL)
+
+                discrete_dist = np.sum((discrete_curve
+                                - np.array([Pnt2d.X(), Pnt2d.Y()]))**2, axis=1)
+
+                proj_idx = np.argmin(discrete_dist)
+
+                discrete_pts_dist[idx] = discrete_dist[proj_idx]
+
+                discrete_pts[idx] = (discrete_curve[proj_idx, 0],
+                                     discrete_curve[proj_idx, 1])
+
+            # Figure out which spline is the closest projection
+            idx = np.argmin(discrete_pts_dist)
+
+            # Save all of the projection information
+            projected_point = gp_Pnt2d(discrete_pts[idx][0],
+                                       discrete_pts[idx][1])
+
+            min_distance = discrete_pts_dist[idx]
+
+            pPnts.append(projected_point)
+            pPara.append(node.parameters[2])
+            pIdx.append(idx)
+            node.corner = False
+
 
         # ==================making sure the pPnts are unique:
         """It happend that somehow the same points were found multiple times"""
@@ -393,6 +450,19 @@ def mesh_by_projecting_nodes_on_BSplineLst(a_BSplineLst, a_nodes, b_BSplineLst,
                     b_nodes[-1].corner = True
                     b_nodes.append(Node(pPnts[0], [LayerID, pIdx[0], pPara[0]]))
 
+            elif len(exterior_corners) == 2 and node.corner == False:
+
+                print("WARNING: Two exterior corners found but node is not seen as a corner.")
+                print("Projection is not included, try increasing crit_angle.")
+                # This case results in one two few b_nodes being added for open
+                # sections and likely a mesh warning later on.
+                # 1. Instead of increasing crit_angle, one could re-evaluate if
+                # 'aglTol' is appropriate.
+                # 2. Or instead one could adopt a different approach here.
+                #
+                # If increasing crit_angle, change the default under
+                # SONATA/cbm/topo/layer/def mesh_layer
+
             # ===CORNERSTYLE 5======
             elif len(exterior_corners) > 2 and node.corner == True:
                 node.cornerstyle = 5
@@ -513,7 +583,6 @@ def mesh_by_projecting_nodes_on_BSplineLst(a_BSplineLst, a_nodes, b_BSplineLst,
 
         except (IndexError):
             print("ERROR:\t IndexError: list index out of range", a_nodes[a])
-            pass
 
         b += 1
 

SONATA/cbm/mesh/mesh_improvements.py

@@ -49,7 +49,12 @@ def modify_sharp_corners(cells, b_BSplineLst, global_minLen, layer_thickness, La
                     FrontNodes = []
                     BackNodes = []
                     distance = (1 + tol) * layer_thickness
-                    for n in MiddleNodes:
+
+                    # Node indices of Middle nodes that fail to project and
+                    # identify front/back.
+                    rm_MiddleNodes = []
+
+                    for mid_idx, n in enumerate(MiddleNodes):
                         pPnts = []
                         pPara = []
                         pIdx = []
@@ -83,11 +88,45 @@ def modify_sharp_corners(cells, b_BSplineLst, global_minLen, layer_thickness, La
                                 BackNodes.append(Node(P, [LayerID, pIdx[i], pPara[i]]))
 
                             else:
-                                print("ERROR: cannot determine FRONT and BACK nodes @ ", c.nodes[0], "because vnp and v01 are orthogonal")
-                                print(pPara)
+                                # This need not be a complete error if other
+                                # Middle nodes have worked.
+                                print("WARNING: cannot determine FRONT and ",
+                                      "BACK nodes @ ",
+                                      c.nodes[0],
+                                      ", so not using this point for sharp ",
+                                      "corner improvement.")
+
+                                rm_MiddleNodes.append(mid_idx)
+
+                        # Make sure 'FrontNodes' and 'BackNodes' are the same
+                        # length and fix.
+                        # This fix is needed in case the front or back projection
+                        # worked, but the other did not. In that case,
+                        # rm_MiddleNodes should have this index.
+                        if mid_idx in rm_MiddleNodes:
+
+                            num_proj = min(len(FrontNodes), len(BackNodes))
+
+                            if len(FrontNodes) - num_proj > 2\
+                                or len(BackNodes) - num_proj > 2:
+                                # Checking every for loop iteration, so these
+                                # lists should not get off by a length of more
+                                # than 1.
+                                print("ERROR: Fixing shaper elements failed @ ", c.nodes[0])
                                 MiddleNodes = []
                                 FrontNodes = []
                                 BackNodes = []
+                            else:
+                                FrontNodes = FrontNodes[:num_proj]
+                                BackNodes = BackNodes[:num_proj]
+
+                    # Remove MiddleNodes that have indices in 'rm_MiddleNodes'
+                    MiddleNodes = [node for k, node in enumerate(MiddleNodes)
+                                   if k not in rm_MiddleNodes]
+
+                    if len(MiddleNodes) == 0 and len(rm_MiddleNodes) > 0:
+                        print("ERROR: All refinements failed for sharp corner @ ",
+                              c.nodes[0], "because vnp and v01 are orthogonal")
 
                     # =====================CREATE FRONT CELLS
                     FrontCellLst = []
@@ -123,8 +162,9 @@ def modify_sharp_corners(cells, b_BSplineLst, global_minLen, layer_thickness, La
                         enhanced_cells.extend(FrontCellLst)
                         enhanced_cells.extend(reversed(BackCellLst))
 
-                        if len(MiddleNodes) == 0:
-                            enhanced_cells.append(c)
+                    if len(MiddleNodes) == 0:
+                        # No refinement found, keep old cell.
+                        enhanced_cells.append(c)
                 else:
                     enhanced_cells.append(c)
 

SONATA/cbm/topo/layer.py

@@ -174,7 +174,9 @@ def determine_a_nodes(self, SegmentLst, global_minLen, display=None):
         self.a_nodes = remove_duplicates_from_list_preserving_order(new_a_nodes)
         self.a_nodes = merge_nodes_if_too_close(self.a_nodes, self.a_BSplineLst, global_minLen, 0.01)
 
-    def mesh_layer(self, SegmentLst, global_minLen, proj_tol_1=9e-2, proj_tol_2=4e-1, crit_angle_1=110, alpha_crit_2=60, growing_factor=1.8, shrinking_factor=0.01, display=None, l0=None):
+    def mesh_layer(self, SegmentLst, global_minLen, proj_tol_1=9e-2, 
+                   proj_tol_2=4e-1, crit_angle_1=150, alpha_crit_2=60, 
+                   growing_factor=1.8, shrinking_factor=0.01, display=None, l0=None):
         """
         The mesh layer function discretizes the layer, which is composed of a 
         a_BsplineLst and a b_BsplineLst. Between the a_BsplineLst and the 
@@ -198,8 +200,13 @@ def mesh_layer(self, SegmentLst, global_minLen, proj_tol_1=9e-2, proj_tol_2=4e-1
         proj_tol_2 = 2e-1:      tolerance value to determine a distance, 
                                 in which the resulting projection point 
                                 has to be. (modify_sharp_corners)
-        crit_angle_1 = 115:     is the critical angle to determine a corner 
-                                if 2 projection points are found.    
+        crit_angle_1 = 150:     is the critical angle to determine a corner 
+                                if 2 projection points are found.
+                                If this value is too small, can hit a case were
+                                both projections are seen as corners, but the
+                                node is not and thus the edge case is not
+                                handled and a mesh warning ultimately gets
+                                raised.
         alpha_crit_2 = 60:      is the critical angle to refine  a corner 
         growing_factor = 1.8:   critical growing factor of cell before 
                                 splitting 

SONATA/cbm/topo/offset.py

@@ -61,8 +61,17 @@ def segment_length(start_idx, end_idx):
             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
+                # Cannot just keep middle point if that would eliminate the
+                # first or last point because that causes the start/end of
+                # the section to move and may cause the segement to no longer
+                # close
+                if i == 0:
+                    combined_points.append(points[0])
                 middle_index = len(close_points) // 2
                 combined_points.append(close_points[middle_index])
+                if j == len(points):
+                    # last point is j-1 by the while loop
+                    combined_points.append(points[-1])
             else:
                 # Keep all points if the length is greater than or equal to the threshold
                 combined_points.extend(close_points)