more fixes

wetb/hawc2/bladeinfo.py → wetb/hawc2/blade.py

@@ -12,6 +12,7 @@ import os
 from wetb.utils.geometry import rad
 from wetb.hawc2.st_file import StFile
 from wetb.hawc2.ae_file import AEFile
+from wetb.hawc2.mainbody import MainBody
 
 
 
@@ -62,7 +63,7 @@ from wetb.hawc2.ae_file import AEFile
 #         """
 
 
-class H2aeroBladeInfo(PCFile, AEFile, AtTimeFile):
+class H2aeroBlade(PCFile, AEFile, AtTimeFile):
     """Provide HAWC2 info about a blade
     
     From AE file:
@@ -302,7 +303,25 @@ class H2aeroBladeInfo(PCFile, AEFile, AtTimeFile):
         return self._Cxxx(radius, alpha, 3, ae_set_nr)
             
             
-class H2BladeInfo(H2aeroBladeInfo):
+    def plot_xz_geometry(self, plt):
+
+        z = np.linspace(self.c2def[0, 2], self.c2def[-1, 2], 100)
+        plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]), label='Center line')
+        plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]) + self.chord(z) / 2, label='Leading edge')
+        plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]) - self.chord(z) / 2, label="Trailing edge")
+        curve_l_nd, x, y, z, twist = self.hawc2_splines()
+        plt.plot(z, x, label='Hawc2spline')
+
+    def plot_yz_geometry(self, plt):
+
+        z = np.linspace(self.c2def[0, 2], self.c2def[-1, 2], 100)
+        plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]), label='Center line')
+        plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) + self.thickness(z) / 100 * self.chord(z) / 2, label='Suction side')
+        plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) - self.thickness(z) / 100 * self.chord(z) / 2, label="Pressure side")
+        curve_l_nd, x, y, z, twist = self.hawc2_splines()
+        plt.plot(z, y, label='Hawc2spline')
+            
+class H2Blade(H2aeroBlade, MainBody):
     """Provide HAWC2 info about a blade
     
     From AE file:
@@ -340,7 +359,10 @@ class H2BladeInfo(H2aeroBladeInfo):
                 blade_name = htcfile.aero.link[2]
             self.mainbody_blade = htcfile.new_htc_structure.get_subsection_by_name(blade_name)
             st_filename = st_filename or os.path.join(htcfile.modelpath, self.mainbody_blade.timoschenko_input.filename[0])
-        H2aeroBladeInfo.__init__(self, htcfile, ae_filename=ae_filename, pc_filename=pc_filename, at_time_filename=at_time_filename, blade_name=blade_name)
+            MainBody.__init__(self, htcfile, blade_name)
+        elif st_filename and os.path.isfile(st_filename):
+            StFile.__init__(self, st_filename)
+        H2aeroBlade.__init__(self, htcfile, ae_filename=ae_filename, pc_filename=pc_filename, at_time_filename=at_time_filename, blade_name=blade_name)
         
 #     def __init__(self, htcfile, ae_filename=None, pc_filename=None, at_time_filename=None, st_filename=None, blade_name=None):
 #         
@@ -353,9 +375,8 @@ class H2BladeInfo(H2aeroBladeInfo):
 #         if os.path.isfile(pc_filename) and os.path.isfile(ae_filename):
 #             PCFile.__init__(self, pc_filename, ae_filename)
 #             blade_radius = self.ae_sets[1][-1,0]
-#         
-        if st_filename and os.path.isfile(st_filename):
-            StFile.__init__(self, st_filename)
+#    
+        
 #         if os.path.isfile(at_time_filename):
 #             AtTimeFile.__init__(self, at_time_filename, blade_radius)
 #             self.curved_length = self.radius_curved_ac()[-1]
@@ -367,6 +388,7 @@ class H2BladeInfo(H2aeroBladeInfo):
 #         self.c2def = np.array([v.values[1:5] for v in mainbody_blade.c2_def if v.name_ == "sec"])
 #         
 #         self.hawc2_splines_data = self.hawc2_splines()       
+    
     @property
     def c2def(self):
         if not hasattr(self, "_c2def"):
@@ -374,7 +396,7 @@ class H2BladeInfo(H2aeroBladeInfo):
         return self._c2def
    
 #         
-# class H2aeroBladeInfo(H2BladeInfo):
+# class H2aeroBlade(H2Blade):
 # 
 #     def __init__(self, at_time_filename, ae_filename, pc_filename, htc_filename):
 #         """

wetb/hawc2/mainbody.py

@@ -12,14 +12,23 @@ from wetb.hawc2.st_file import StFile
 
 
 class MainBody():
-    def __init__(self, htc_filename, modelpath, body_name):
-        self.htcfile = htcfile = HTCFile(htc_filename, modelpath)
+    def __init__(self, htcfile, body_name):
+        if isinstance(htcfile, str):
+            htcfile = HTCFile(htcfile)
+        self.htcfile = htcfile
         s = htcfile.new_htc_structure
         main_bodies = {s[k].name[0]:s[k] for k in s.keys() if s[k].name_ == "main_body"}
         self.main_body = main_bodies[body_name]
+        self.main_body = s.get_subsection_by_name(body_name)
         self.stFile = StFile(os.path.join(htcfile.modelpath, self.main_body.timoschenko_input.filename[0]))
-        self.c2def = np.array([v.values[1:5] for v in self.main_body.c2_def if v.name_ == "sec"])
+        #self.c2def = np.array([v.values[1:5] for v in self.main_body.c2_def if v.name_ == "sec"])
         self.concentrated_mass = [cm.values for cm in self.main_body if cm.name_.startswith('concentrated_mass')]
+        
+    @property
+    def c2def(self):
+        if not hasattr(self, "_c2def"):
+            self._c2def = np.array([v.values[1:5] for v in self.main_body.c2_def if v.name_ == "sec"])
+        return self._c2def
 
     def plot_xz_geometry(self, plt=None):
         if plt is None:
@@ -43,7 +52,7 @@ class MainBody():
         plt.xlabel("z")
         plt.ylabel("y")
         z = np.linspace(self.c2def[0, 2], self.c2def[-1, 2], 100)
-        plt.plot(self.c2def[:, 2], self.c2def[:, 1], label='Center line')
+        plt.plot(self.c2def[:, 2], self.c2def[:, 1], '.-', label='Center line')
         plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) + self.stFile.y_e(z), label='Elastic center')
         plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) + self.stFile.y_cg(z), label='Mass center')
         plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) + self.stFile.y_sh(z), label='Shear center')
@@ -52,106 +61,106 @@ class MainBody():
         plt.legend()
 
         
-
-
-class BladeData(object):
-    def plot_xz_geometry(self, plt):
-
-        z = np.linspace(self.c2def[0, 2], self.c2def[-1, 2], 100)
-        plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]), label='Center line')
-        plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]) + self.pcFile.chord(z) / 2, label='Leading edge')
-        plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]) - self.pcFile.chord(z) / 2, label="Trailing edge")
-        x, y, z = self.hawc2_splines()
-        #plt.plot(z, x, label='Hawc2spline')
-
-    def plot_yz_geometry(self, plt):
-
-        z = np.linspace(self.c2def[0, 2], self.c2def[-1, 2], 100)
-        plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]), label='Center line')
-        plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) + self.pcFile.thickness(z) / 100 * self.pcFile.chord(z) / 2, label='Suction side')
-        plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) - self.pcFile.thickness(z) / 100 * self.pcFile.chord(z) / 2, label="Pressure side")
-        x, y, z = self.hawc2_splines()
-        #plt.plot(z, y, label='Hawc2spline')
-
-    def hawc2_splines(self):
-        curve_z = np.r_[0, np.cumsum(np.sqrt(np.sum(np.diff(self.c2def[:, :3], 1, 0) ** 2, 1)))]
-        curve_z_nd = curve_z / curve_z[-1]
-
-        def akima(x, y):
-            n = len(x)
-            var = np.zeros((n + 3))
-            z = np.zeros((n))
-            co = np.zeros((n, 4))
-            for i in range(n - 1):
-                var[i + 2] = (y[i + 1] - y[i]) / (x[i + 1] - x[i])
-            var[n + 1] = 2 * var[n] - var[n - 1]
-            var[n + 2] = 2 * var[n + 1] - var[n]
-            var[1] = 2 * var[2] - var[3]
-            var[0] = 2 * var[1] - var[2]
-
-            for i in range(n):
-                wi1 = abs(var[i + 3] - var[i + 2])
-                wi = abs(var[i + 1] - var[i])
-                if (wi1 + wi) == 0:
-                    z[i] = (var[i + 2] + var[i + 1]) / 2
-                else:
-                    z[i] = (wi1 * var[i + 1] + wi * var[i + 2]) / (wi1 + wi)
-
-            for i in range(n - 1):
-                dx = x[i + 1] - x[i]
-                a = (z[i + 1] - z[i]) * dx
-                b = y[i + 1] - y[i] - z[i] * dx
-                co[i, 0] = y[i]
-                co[i, 1] = z[i]
-                co[i, 2] = (3 * var[i + 2] - 2 * z[i] - z[i + 1]) / dx
-                co[i, 3] = (z[i] + z[i + 1] - 2 * var[i + 2]) / dx ** 2
-            co[n - 1, 0] = y[n - 1]
-            co[n - 1, 1] = z[n - 1]
-            co[n - 1, 2] = 0
-            co[n - 1, 3] = 0
-            return co
-
-        def coef2spline(s, co):
-            x, y = [], []
-            for i, c in enumerate(co.tolist()[:-1]):
-                p = np.poly1d(c[::-1])
-                z = np.linspace(0, s[i + 1] - s[i ], 10)
-                x.extend(s[i] + z)
-                y.extend(p(z))
-            return y
-
-        x, y, z = [coef2spline(curve_z_nd, akima(curve_z_nd, self.c2def[:, i])) for i in range(3)]
-        return x, y, z
-
-
-
-class Blade(MainBody, BladeData):
-    def __init__(self, htc_filename, modelpath=None, blade_number=1):
-        
-        self.htcfile = htcfile = HTCFile(htc_filename, modelpath)
-        
-        blade_name = [link[2] for link in htcfile.aero if link.name_.startswith('link') and link[0]==blade_number][0]
-        MainBody.__init__(self, htc_filename, modelpath, blade_name)
-        self.pcFile = PCFile(os.path.join(htcfile.modelpath, htcfile.aero.pc_filename[0]),
-                        os.path.join(htcfile.modelpath, htcfile.aero.ae_filename[0]))
-        
-    def plot_xz_geometry(self, plt=None):
-        if plt is None:
-            import matplotlib.pyplot as plt
-            plt.figure()
-
-        MainBody.plot_xz_geometry(self, plt)
-        BladeData.plot_xz_geometry(self, plt=plt)
-        plt.legend()
-    
-    def plot_geometry_yz(self, plt=None):
-        if plt is None:
-            import matplotlib.pyplot as plt
-            plt.figure()
-
-        BladeData.plot_yz_geometry(self, plt=plt)
-        MainBody.plot_yz_geometry(self, plt)
-        
+# 
+# 
+# class BladeData(object):
+#     def plot_xz_geometry(self, plt):
+# 
+#         z = np.linspace(self.c2def[0, 2], self.c2def[-1, 2], 100)
+#         plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]), label='Center line')
+#         plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]) + self.pcFile.chord(z) / 2, label='Leading edge')
+#         plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]) - self.pcFile.chord(z) / 2, label="Trailing edge")
+#         x, y, z = self.hawc2_splines()
+#         #plt.plot(z, x, label='Hawc2spline')
+# 
+#     def plot_yz_geometry(self, plt):
+# 
+#         z = np.linspace(self.c2def[0, 2], self.c2def[-1, 2], 100)
+#         plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]), label='Center line')
+#         plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) + self.pcFile.thickness(z) / 100 * self.pcFile.chord(z) / 2, label='Suction side')
+#         plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) - self.pcFile.thickness(z) / 100 * self.pcFile.chord(z) / 2, label="Pressure side")
+#         x, y, z = self.hawc2_splines()
+#         #plt.plot(z, y, label='Hawc2spline')
+# 
+#     def hawc2_splines(self):
+#         curve_z = np.r_[0, np.cumsum(np.sqrt(np.sum(np.diff(self.c2def[:, :3], 1, 0) ** 2, 1)))]
+#         curve_z_nd = curve_z / curve_z[-1]
+# 
+#         def akima(x, y):
+#             n = len(x)
+#             var = np.zeros((n + 3))
+#             z = np.zeros((n))
+#             co = np.zeros((n, 4))
+#             for i in range(n - 1):
+#                 var[i + 2] = (y[i + 1] - y[i]) / (x[i + 1] - x[i])
+#             var[n + 1] = 2 * var[n] - var[n - 1]
+#             var[n + 2] = 2 * var[n + 1] - var[n]
+#             var[1] = 2 * var[2] - var[3]
+#             var[0] = 2 * var[1] - var[2]
+# 
+#             for i in range(n):
+#                 wi1 = abs(var[i + 3] - var[i + 2])
+#                 wi = abs(var[i + 1] - var[i])
+#                 if (wi1 + wi) == 0:
+#                     z[i] = (var[i + 2] + var[i + 1]) / 2
+#                 else:
+#                     z[i] = (wi1 * var[i + 1] + wi * var[i + 2]) / (wi1 + wi)
+# 
+#             for i in range(n - 1):
+#                 dx = x[i + 1] - x[i]
+#                 a = (z[i + 1] - z[i]) * dx
+#                 b = y[i + 1] - y[i] - z[i] * dx
+#                 co[i, 0] = y[i]
+#                 co[i, 1] = z[i]
+#                 co[i, 2] = (3 * var[i + 2] - 2 * z[i] - z[i + 1]) / dx
+#                 co[i, 3] = (z[i] + z[i + 1] - 2 * var[i + 2]) / dx ** 2
+#             co[n - 1, 0] = y[n - 1]
+#             co[n - 1, 1] = z[n - 1]
+#             co[n - 1, 2] = 0
+#             co[n - 1, 3] = 0
+#             return co
+# 
+#         def coef2spline(s, co):
+#             x, y = [], []
+#             for i, c in enumerate(co.tolist()[:-1]):
+#                 p = np.poly1d(c[::-1])
+#                 z = np.linspace(0, s[i + 1] - s[i ], 10)
+#                 x.extend(s[i] + z)
+#                 y.extend(p(z))
+#             return y
+# 
+#         x, y, z = [coef2spline(curve_z_nd, akima(curve_z_nd, self.c2def[:, i])) for i in range(3)]
+#         return x, y, z
+
+
+
+# class Blade(MainBody, BladeData):
+#     def __init__(self, htc_filename, modelpath=None, blade_number=1):
+#         
+#         self.htcfile = htcfile = HTCFile(htc_filename, modelpath)
+#         
+#         blade_name = [link[2] for link in htcfile.aero if link.name_.startswith('link') and link[0]==blade_number][0]
+#         MainBody.__init__(self, htc_filename, modelpath, blade_name)
+#         self.pcFile = PCFile(os.path.join(htcfile.modelpath, htcfile.aero.pc_filename[0]),
+#                         os.path.join(htcfile.modelpath, htcfile.aero.ae_filename[0]))
+#         
+#     def plot_xz_geometry(self, plt=None):
+#         if plt is None:
+#             import matplotlib.pyplot as plt
+#             plt.figure()
+# 
+#         MainBody.plot_xz_geometry(self, plt)
+#         BladeData.plot_xz_geometry(self, plt=plt)
+#         plt.legend()
+#     
+#     def plot_geometry_yz(self, plt=None):
+#         if plt is None:
+#             import matplotlib.pyplot as plt
+#             plt.figure()
+# 
+#         BladeData.plot_yz_geometry(self, plt=plt)
+#         MainBody.plot_yz_geometry(self, plt)
+#         
 
 
 

wetb/hawc2/tests/test_blade.py

+'''
+Created on 3. maj 2017
+
+@author: mmpe
+'''
+import os
+import unittest
+
+from wetb.hawc2.blade import H2Blade
+import numpy as np
+
+tfp = os.path.join(os.path.dirname(__file__), 'test_files/')  # test file path
+class Test(unittest.TestCase):
+
+
+    def testBladeInfo(self):
+        bi = H2Blade(tfp + "simulation_setup/DTU10MWRef6.0/htc/DTU_10MW_RWT.htc")
+        if 0:
+            import matplotlib.pyplot as plt
+            print (dir(bi))
+            #print (bi.radius_s())
+            plt.plot(bi.radius_s(), bi.twist())
+            plt.plot(bi.c2def[:,2], bi.c2def[:,3])
+            x = np.linspace(0,1,1000)
+            plt.plot(bi.blade_radius*x, bi.c2nd(x)[:,3])
+            plt.show()
+
+            
+    def testBladeInfo_AE(self):
+        bi = H2Blade(None, tfp + "NREL_5MW_ae.txt", tfp + "NREL_5MW_pc.txt")
+        self.assertEqual(bi.thickness(32), 23.78048780487805)
+        self.assertEqual(bi.chord(32), 3.673)
+        self.assertEqual(bi.pc_set_nr(32), 1)
+
+    def testBladeInfo_PC_AE(self):
+        bi = H2Blade(None, tfp + "NREL_5MW_ae.txt", tfp + "NREL_5MW_pc.txt")
+        self.assertEqual(bi.CL(36, 10), 1.358)
+        self.assertEqual(bi.CD(36, 10), 0.0255)
+        self.assertEqual(bi.CM(36, 10), -0.1103)
+            
+#     def test_curved_length2radius(self):
+#         bi = H2BladeInfo(tfp + "simulation_setup/DTU10MWRef6.0/htc/DTU_10MW_RWT.htc")
+#         np.testing.assert_array_equal(bi.xyztwist(86.4979, curved_length=True), bi.xyztwist(86.3655))
+        
+            
+if __name__ == "__main__":
+    #import sys;sys.argv = ['', 'Test.testBladeInfo']
+    unittest.main()
\ No newline at end of file

wetb/hawc2/tests/test_main_body.py

@@ -5,13 +5,14 @@ Created on 01/08/2016
 '''
 import unittest
 import os
-from wetb.hawc2.mainbody import MainBody, Blade
+from wetb.hawc2.mainbody import MainBody
+from wetb.hawc2.blade import H2Blade
 
 tfp = os.path.join(os.path.dirname(__file__), 'test_files/')  # test file path
 class TestMainBody(unittest.TestCase):
     
     def test_MainBody(self):
-        mainbody = MainBody(tfp+"htcfiles/test.htc", "../", "towertop")
+        mainbody = MainBody(tfp+"htcfiles/test.htc", "blade1")
         if 0:
             import matplotlib.pyplot as plt
             plt.figure()
@@ -21,7 +22,7 @@ class TestMainBody(unittest.TestCase):
             plt.show()
              
     def test_Blade(self):
-        blade = Blade(tfp+"htcfiles/test.htc", "../")
+        blade = H2Blade(tfp+"htcfiles/test.htc")
         if 0:
             import matplotlib.pyplot as plt
             plt.figure()