From 7acb49a00b5fbaf0e926f86d5bb92c7a892e2368 Mon Sep 17 00:00:00 2001
From: Malo Rosemeier <malo.rosemeier@iwes.fraunhofer.de>
Date: Fri, 21 Jul 2017 11:46:20 +0200
Subject: [PATCH] renaming and clean up

---
 ...ip_blade.py => finstrip_bladestability.py} |  77 ++++++----
 ...ade.py => test_finstrip_bladestability.py} | 135 +++++++++++++-----
 2 files changed, 152 insertions(+), 60 deletions(-)
 rename finstrip_wrapper/omdao/{finstrip_blade.py => finstrip_bladestability.py} (92%)
 rename finstrip_wrapper/omdao/test/{test_finstrip_blade.py => test_finstrip_bladestability.py} (66%)

diff --git a/finstrip_wrapper/omdao/finstrip_blade.py b/finstrip_wrapper/omdao/finstrip_bladestability.py
similarity index 92%
rename from finstrip_wrapper/omdao/finstrip_blade.py
rename to finstrip_wrapper/omdao/finstrip_bladestability.py
index ce2c380..852f766 100644
--- a/finstrip_wrapper/omdao/finstrip_blade.py
+++ b/finstrip_wrapper/omdao/finstrip_bladestability.py
@@ -6,14 +6,11 @@ from PGL.components.airfoil import AirfoilShape
 from becas_wrapper.becas_bladestructure import set_cs_size
 
 
-class FINSTRIPCSStructure(Component):
+class FINSTRIPCSStability(Component):
     """
-    Component for computing beam structural properties
-    using the cross-sectional structure code FINSTRIP.
+    Component for computing blade stability resistance using FINSTRIP.
 
-    The code firstly calls CS2DtoFINSTRIP which is a wrapper around
-    shellexpander that comes with FINSTRIP, and second
-    calls FINSTRIP using a file interface.
+    The code calls FINSTRIP using a file interface.
 
     parameters
     ----------
@@ -52,7 +49,7 @@ class FINSTRIPCSStructure(Component):
         ni_chord: int
             number of points defining the cross-section shape
         """
-        super(FINSTRIPCSStructure, self).__init__()
+        super(FINSTRIPCSStability, self).__init__()
 
         self.basedir = os.getcwd()
         self.finstrip_hash = finstrip_hash
@@ -269,11 +266,14 @@ class FINSTRIPCSStructure(Component):
         os.chdir(self.basedir)
 
 
-class FINSTRIPCSStructureLoader(FINSTRIPCSStructure):
+class FINSTRIPCSStabilityLoader(FINSTRIPCSStability):
+    '''
+    NOTE: not yet tested
+    '''
 
     def __init__(self, name, finstrip_hash, config, st3d, s, ni_chord, cs_size, cs_size_ref, pickled_unkn,
                  grp_name=''):
-        super(FINSTRIPCSStructureLoader, self).__init__(
+        super(FINSTRIPCSStabilityLoader, self).__init__(
             name, finstrip_hash, config, st3d, s, ni_chord, cs_size, cs_size_ref)
 
         self.pickled_unkn = pickled_unkn
@@ -302,23 +302,35 @@ class Slice(Component):
     blade_surface_norm_st: array
         blade surface with structural discretization, no twist and prebend.
         Size: ((ni_chord, nsec, 3))
+    blade_length: float
+        blade length
 
     returns
     -------
-    sec<xxx>DPs: array
+    stab:sec<xxx>DPs: array
         Vector of DPs along chord for each section. Size (nDP)
-    sec<xxx>coords: array
+    stab:sec<xxx>coords: array
         Array of cross section coords shapes. Size ((ni_chord, 3))
+    stab:sec<xxx>tvec: array
+        Array of region thicknesses and angles. Size (len(self._varnames) * 2)        
+    stab:sec<xxx>rot_z_st: float
+        geometric twist of cross section
+    stab:sec<xxx>csprops: array (optional)
+        Array of cross section props as of BECAS. Size ((cs_size))        
+    stab:sec<xxx>csprops_ref: array (optional)
+        Array of cross section props ref as of BECAS. Size ((cs_size_ref))
     """
 
     def __init__(self, st3d, sdim, config):
         """
         parameters
         ----------
-        DPs: array
-            DPs array, size: ((nsec, nDP))
+        st3d: dict
+            dictionary of blade structure properties
         sdim: array
             blade surface. Size: ((ni_chord, nsec, 3))
+        config: dict
+            configuration
         """
         super(Slice, self).__init__()
 
@@ -418,6 +430,9 @@ class Slice(Component):
 
 
 class SliceLoader(Slice):
+    '''
+    NOTE: not yet tested
+    '''
 
     def __init__(self, st3d, sdim, pickled_unkn, grp_name=''):
         super(SliceLoader, self).__init__(st3d, sdim)
@@ -444,12 +459,12 @@ class PostprocessCS(Component):
     parameters
     ----------
     f_eig<xxx>: array
-        array of cross section eigen values. size (nmode,2)
+        array of cross section eigen values. size (ncase,nmode,2)
 
     returns
     -------
     blade_f_eig: array
-        array of buckling resistances. Size ((nsec,nmode,2)).
+        array of buckling resistances. Size ((ncase,nsec,nmode,2)).
     """
 
     def __init__(self, nsec, nmode, ncase):
@@ -458,6 +473,10 @@ class PostprocessCS(Component):
         ----------
         nsec: int
             number of blade sections.
+        nmode: int
+            number of buckling modes.
+        ncase: int
+            number of load cases.
         """
         super(PostprocessCS, self).__init__()
 
@@ -479,6 +498,9 @@ class PostprocessCS(Component):
 
 
 class PostprocessCSLoader(PostprocessCS):
+    '''
+    NOTE: not yet tested
+    '''
 
     def __init__(self, nsec, cs_size, cs_size_ref, pickled_unkn, grp_name='', promoted=[]):
         super(PostprocessCSLoader, self).__init__(nsec, cs_size, cs_size_ref)
@@ -497,7 +519,7 @@ class PostprocessCSLoader(PostprocessCS):
                     unknowns[ku] = vp
 
 
-class FINSTRIPBeamStructure(Group):
+class FINSTRIPBladeStability(Group):
     """
     Group for computing blade buckling resistance using FINSTRIP.
 
@@ -527,7 +549,7 @@ class FINSTRIPBeamStructure(Group):
     returns
     -------
     blade_f_eig: array
-        array of buckling resistances. Size ((nsec,nmode,2)).
+        array of buckling resistances. Size ((ncase,nsec,nmode,2)).
     """
 
     def __init__(self, group, config, st3d, sdim, ncase):
@@ -537,6 +559,8 @@ class FINSTRIPBeamStructure(Group):
 
         parameters
         ----------
+        group: omdao Group
+            group FINSTRIPBladeStability is inside
         config: dict
             dictionary of inputs for the cs_code class
         st3d: dict
@@ -544,8 +568,10 @@ class FINSTRIPBeamStructure(Group):
         surface: array
             blade surface with structural discretization.
             Size: ((ni_chord, nsec, 3))
+        ncase: int
+            number of load cases.
         """
-        super(FINSTRIPBeamStructure, self).__init__()
+        super(FINSTRIPBladeStability, self).__init__()
 
         # check that the config is ok
         if not 'FINSTRIPWrapper' in config.keys():
@@ -582,7 +608,7 @@ class FINSTRIPBeamStructure(Group):
 
         for i in range(nsec):
             secname = 'sec%03d' % i
-            par.add(secname, FINSTRIPCSStructure(secname, self.__hash__(), config, st3d,
+            par.add(secname, FINSTRIPCSStability(secname, self.__hash__(), config, st3d,
                                                  st3d['s'][i], sdim[0], ncase), promotes=['*'])
 
         promotions = ['blade_f_eig']
@@ -593,17 +619,20 @@ class FINSTRIPBeamStructure(Group):
             self.connect('%s:f_eig' % secname, 'postpro.f_eig%03d' % i)
 
 
-class FINSTRIPBeamStructureLoader(Group):
+class FINSTRIPBladeStabilityLoader(Group):
     """
     Component for loading eigen values stored in a database by a
     previous run.
 
-    This component can be replaced by FINSTRIPBeamStructure group in case the FINSTRIP
+    This component can be replaced by FINSTRIPBladeStability group in case the FINSTRIP
     run shall be skipped, i.e. when the structure or load case does not change during 
     iterations.
 
     Note: If parameters are promoted on the level higher than the group 
-    FINSTRIPBeamStructure belongs to, they need to be listed in the promoted list.
+    FINSTRIPBladeStability belongs to, they need to be listed in the promoted list.
+
+    NOTE: not yet tested
+
     """
 
     def __init__(self, group, config, st3d, sdim, pickled_unkn, grp_name=''):
@@ -631,7 +660,7 @@ class FINSTRIPBeamStructureLoader(Group):
             >>> pickled_unkn = data['Unknowns']
 
         """
-        super(FINSTRIPBeamStructureLoader, self).__init__()
+        super(FINSTRIPBladeStabilityLoader, self).__init__()
 
         # check that the config is ok
         if not 'FINSTRIPWrapper' in config.keys():
@@ -668,7 +697,7 @@ class FINSTRIPBeamStructureLoader(Group):
 
         for i in range(nsec):
             secname = 'sec%03d' % i
-            par.add(secname, FINSTRIPCSStructureLoader(secname, self.__hash__(), config, st3d,
+            par.add(secname, FINSTRIPCSStabilityLoader(secname, self.__hash__(), config, st3d,
                                                        st3d['s'][i], sdim[0], pickled_unkn, grp_name), promotes=['*'])
 
         promotions = ['blade_f_eig']
diff --git a/finstrip_wrapper/omdao/test/test_finstrip_blade.py b/finstrip_wrapper/omdao/test/test_finstrip_bladestability.py
similarity index 66%
rename from finstrip_wrapper/omdao/test/test_finstrip_blade.py
rename to finstrip_wrapper/omdao/test/test_finstrip_bladestability.py
index f01b045..e6a7371 100644
--- a/finstrip_wrapper/omdao/test/test_finstrip_blade.py
+++ b/finstrip_wrapper/omdao/test/test_finstrip_bladestability.py
@@ -18,15 +18,12 @@ from fusedwind.turbine.structure import read_bladestructure, \
     SplinedBladeStructure, \
     SplinedBladeStructureParam2
 
-from becas_wrapper.becas_bladestructure import BECASBeamStructureKM,\
-    BECASBeamStructure
-from becas_wrapper.becas_stressrecovery import BECASStressRecovery
+from becas_wrapper.becas_bladestructure import BECASBeamStructure
 
-from distutils.spawn import find_executable
-from finstrip_wrapper.omdao.finstrip_blade import FINSTRIPBeamStructure
+from finstrip_wrapper.omdao.finstrip_bladestability import FINSTRIPBladeStability
 from robe.lib.load import convert_point_forces_to_moment_distro
 
-_matlab_installed = find_executable('matlab')
+_becas_installed = os.getenv('BECAS_BASEDIR')
 
 if MPI:
     # if you called this script with 'mpirun', then use the petsc data passing
@@ -36,7 +33,7 @@ else:
     from openmdao.core.basic_impl import BasicImpl as impl
 
 try:
-    DRYRUN = os.environ['BECAS_DRYRUN']
+    DRYRUN = os.environ['FINSTRIP_DRYRUN']
 except:
     DRYRUN = '0'
 if DRYRUN == '1':
@@ -44,8 +41,10 @@ if DRYRUN == '1':
 else:
     DRYRUN = False
 
+PLOT_MODE_SHAPES = True
 
-def configure_BECASBeamStructure(nsec, exec_mode, path_data, dry_run=False, FPM=False, with_sr=False, param2=False, with_becas=False):
+
+def configure_FINSTRIPBladeStability(nsec, exec_mode, path_data, dry_run=False, FPM=False, with_sr=False, param2=False, with_becas=False):
 
     p = Problem(impl=impl, root=Group())
 
@@ -126,7 +125,7 @@ def configure_BECASBeamStructure(nsec, exec_mode, path_data, dry_run=False, FPM=
         #        config, s_st, 2), promotes=['*'])
 
     cfg = {}
-    cfg['nmode'] = 20
+    cfg['nmode'] = 20  # number of buckling modes
     cfg['core_mats'] = ['balsa']
     cfg['elsize_factor'] = 0.015
     cfg['lambda_min'] = 0.01  # minimum half wavelength in m
@@ -135,8 +134,8 @@ def configure_BECASBeamStructure(nsec, exec_mode, path_data, dry_run=False, FPM=
     config['FINSTRIPWrapper'] = cfg
 
     ncase = 2
-    p.root.add('stability', FINSTRIPBeamStructure(p.root, config, st3dn,
-                                                  (200, nsec_st, 3), ncase), promotes=['*'])
+    p.root.add('stability', FINSTRIPBladeStability(p.root, config, st3dn,
+                                                   (200, nsec_st, 3), ncase), promotes=['*'])
 
     p.setup()
     for k, v in pf.iteritems():
@@ -190,6 +189,61 @@ def configure_BECASBeamStructure(nsec, exec_mode, path_data, dry_run=False, FPM=
     return p
 
 
+def plot_mode_shapes(p, figname):
+    ncase = p['load_cases_sec%03d' % 0].shape[0]
+    nsec = len(p['s_st'])
+    f_eig_min = np.zeros((ncase, nsec))
+    f_eig_min_idx = np.zeros((ncase, nsec)).astype(int)
+    for case in range(ncase):
+        for sec in range(nsec):
+            f_eig_min[case, sec] = np.min(
+                p['blade_f_eig'][case, sec, :, 1])
+            f_eig_min_idx[case, sec] = p[
+                'blade_f_eig'][case, sec, :, 1].argmin()
+
+    import matplotlib.pylab as plt
+    fig, ax = plt.subplots()
+    sf = 500
+    case_cols = ['r', 'g']
+    for case in range(ncase):
+        for sec in range(nsec):
+            # plot nodes
+            nnode_surf = p['sec%03d:nnode_surf' % sec]
+            node_loc_surf = p['sec%03d:node_loc_surf' %
+                              sec][:nnode_surf, :]
+            ax.plot(node_loc_surf[:, 0], node_loc_surf[:, 1], 'b-')
+
+            nnode_web = p['sec%03d:nnode_web' % sec]
+            nweb = len(nnode_web)
+            node_loc_web = p['sec%03d:node_loc_web' %
+                             sec][:nweb, :, :]
+            for w in range(nweb):
+                plt.plot(
+                    node_loc_web[w, :nnode_web[w], 0],
+                    node_loc_web[w, :nnode_web[w], 1], 'b-')
+
+            # plot displacements
+            node_u_surf = p['sec%03d:node_u_surf' %
+                            sec][case, f_eig_min_idx[case, sec], :nnode_surf, :]
+            node_u_web = p['sec%03d:node_u_web' %
+                           sec][case, :, f_eig_min_idx[case, sec], :, :]
+
+            ax.plot(node_loc_surf[:, 0] + sf * node_u_surf[:, 0],
+                    node_loc_surf[:, 1] + sf * node_u_surf[:, 1],
+                    case_cols[case] + '-')
+
+            for w in range(nweb):
+                ax.plot(node_loc_web[w, :nnode_web[w], 0] +
+                        sf * node_u_web[w, :nnode_web[w], 0],
+                        node_loc_web[w, :nnode_web[w], 1] +
+                        sf * node_u_web[w, :nnode_web[w], 1],
+                        case_cols[case] + '-')
+
+    plt.axis('equal')
+    fig.savefig(os.path.join(figname + '_mode_shapes' + '.pdf'))
+    # plt.show()
+
+
 class FINSTRIPBeamStructureTestCase(unittest.TestCase):
 
     def setUp(self):
@@ -199,44 +253,50 @@ class FINSTRIPBeamStructureTestCase(unittest.TestCase):
         pass
 
     #@unittest.skip("skipping")
-    @unittest.skipIf(not _matlab_installed,
-                     "Matlab not available on this system")
-    def test_standard_matlab(self):
-
-        p = configure_BECASBeamStructure(4, 'matlab', 'data',
-                                         dry_run=DRYRUN,
-                                         FPM=False,
-                                         with_sr=True,
-                                         param2=False,
-                                         with_becas=True)
-        p.run()
+    def test_finstrip(self):
 
-        ncase = 2
+        p = configure_FINSTRIPBladeStability(4, 'octave', 'data',
+                                             dry_run=DRYRUN,
+                                             FPM=False,
+                                             with_sr=True,
+                                             param2=False,
+                                             with_becas=False)
+
+        p.run()
+        ncase = p['load_cases_sec%03d' % 0].shape[0]
         nsec = len(p['s_st'])
         f_eig_min = np.zeros((ncase, nsec))
+        f_eig_min_idx = np.zeros((ncase, nsec)).astype(int)
         for case in range(ncase):
             for sec in range(nsec):
                 f_eig_min[case, sec] = np.min(
                     p['blade_f_eig'][case, sec, :, 1])
+                f_eig_min_idx[case, sec] = p[
+                    'blade_f_eig'][case, sec, :, 1].argmin()
 
-        f_eig_min_exp = np.array([[0.727,  0.833,  1.385,  1.],
-                                  [4.734,  3.066,  1.694,  1.]])
+        f_eig_min_exp = np.array([[0.996,  1.877,  2.09,  1.],
+                                  [6.484,  5.243,  2.011,  1.]])
         for case in range(ncase):
             for i, f in enumerate(f_eig_min[case, :]):
                 self.assertAlmostEqual(
                     f, f_eig_min_exp[case, i], places=3)
 
-    #@unittest.skip("skipping")
-    def test_finstrip(self):
-
-        p = configure_BECASBeamStructure(4, 'matlab', 'data',
-                                         dry_run=DRYRUN,
-                                         FPM=False,
-                                         with_sr=True,
-                                         param2=False,
-                                         with_becas=False)
+        if PLOT_MODE_SHAPES:
+            plot_mode_shapes(p, 'finstrip')
 
+    #@unittest.skip("skipping")
+    @unittest.skipIf(not _becas_installed,
+                     "BECAS not available on this system")
+    def test_finstrip_becas(self):
+
+        p = configure_FINSTRIPBladeStability(4, 'octave', 'data',
+                                             dry_run=DRYRUN,
+                                             FPM=False,
+                                             with_sr=True,
+                                             param2=False,
+                                             with_becas=True)
         p.run()
+
         ncase = 2
         nsec = len(p['s_st'])
         f_eig_min = np.zeros((ncase, nsec))
@@ -245,13 +305,16 @@ class FINSTRIPBeamStructureTestCase(unittest.TestCase):
                 f_eig_min[case, sec] = np.min(
                     p['blade_f_eig'][case, sec, :, 1])
 
-        f_eig_min_exp = np.array([[0.996,  1.877,  2.09,  1.],
-                                  [6.484,  5.243,  2.011,  1.]])
+        f_eig_min_exp = np.array([[0.727,  0.833,  1.385,  1.],
+                                  [4.734,  3.066,  1.694,  1.]])
         for case in range(ncase):
             for i, f in enumerate(f_eig_min[case, :]):
                 self.assertAlmostEqual(
                     f, f_eig_min_exp[case, i], places=3)
 
+        if PLOT_MODE_SHAPES:
+            plot_mode_shapes(p, 'finstrip_becas')
+
 
 if __name__ == "__main__":
 
-- 
GitLab