Added bladeData.py + test

wetb/hawc2/bladeData.py

+'''
+Created on 01/08/2016
+
+@author: MMPE
+'''
+from wetb.hawc2.pc_file import PCFile
+from wetb.hawc2.ae_file import AEFile
+from wetb.hawc2.htc_file import HTCFile
+import os
+import numpy as np
+from wetb.hawc2.st_file import StFile
+
+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 H2BladeData(BladeData):
+    def __init__(self, htc_filename, modelpath):
+        self.htcfile = htcfile = HTCFile(htc_filename, modelpath)
+        self.pcFile = PCFile(os.path.join(htcfile.modelpath, htcfile.aero.pc_filename[0]),
+                        os.path.join(htcfile.modelpath, htcfile.aero.ae_filename[0]))
+        s = htcfile.new_htc_structure
+        blade_name = htcfile.aero.link[2]
+        mainbodies = [s[k] for k in s.keys() if s[k].name_ == "main_body"]
+
+        blade_main_body = [mb for mb in mainbodies if mb.name[0] == blade_name][0]
+        self.stFile = StFile(os.path.join(htcfile.modelpath, blade_main_body.timoschenko_input.filename[0]))
+        self.c2def = np.array([v.values[1:5] for v in blade_main_body.c2_def if v.name_ == "sec"])
+
+    def plot_geometry(self, plt=None):
+        if plt is None:
+            import matplotlib.pyplot as plt
+
+        BladeData.plot_xz_geometry(self, plt=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]) + self.stFile.x_e(z), label='Elastic center')
+        plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]) + self.stFile.x_cg(z), label='Mass center')
+        plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]) + self.stFile.x_sh(z), label='Shear center')
+
+        plt.legend()
+        plt.figure()
+        BladeData.plot_yz_geometry(self, plt=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]) + 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')
+        plt.legend()
+
+        plt.show()
+
+
+
+
+class H2AeroBladeData(BladeData):
+    def __init__(self, htc_filename, modelpath):
+        self.htcfile = htcfile = HTCFile(htc_filename, modelpath)
+        self.pcFile = PCFile(os.path.join(htcfile.modelpath, htcfile.aero.pc_filename[0]),
+                        os.path.join(htcfile.modelpath, htcfile.aero.ae_filename[0]))
+        self.c2def = np.array([v.values[1:5] for v in htcfile.blade_c2_def if v.name_ == "sec"])
+
+    def plot_geometry(self, plt=None):
+        if plt is None:
+            import matplotlib.pyplot as plt
+
+        BladeData.plot_xz_geometry(self, plt=plt)
+        z = np.linspace(self.c2def[0, 2], self.c2def[-1, 2], 100)
+        plt.legend()
+        plt.figure()
+        BladeData.plot_yz_geometry(self, plt=plt)
+        z = np.linspace(self.c2def[0, 2], self.c2def[-1, 2], 100)
+        plt.legend()
+
+        plt.show()
+
+

wetb/hawc2/tests/test_bladeData.py

+'''
+Created on 01/08/2016
+
+@author: MMPE
+'''
+import unittest
+import os
+from wetb.hawc2.bladeData import H2AeroBladeData, H2BladeData
+tfp = os.path.join(os.path.dirname(__file__), 'test_files/')  # test file path
+class TestBladeData(unittest.TestCase):
+    
+    def testBladedata(self):
+        #bd = H2AeroBladeData(tfp + "h2aero_tb/htc/h2aero.htc", '../')
+        #bd = H2AeroBladeData(r'C:\mmpe\programming\python\Phd\pitot_tube\tests\test_files\swt36_107\h2a\htc\templates/swt3.6_107.htc', "../../")
+        #print (bd.pcFile.chord(36))
+        #bd = H2BladeData(r"C:\mmpe\HAWC2\models\SWT3.6-107\original_data\SWT4.0-130\dlc12_wsp06_wdir010_s18002.htc", ".")
+        bd = H2BladeData(r"C:\mmpe\HAWC2\models\SWT3.6-107\htc/stat6/stat6_10.0_0.htc", "../../")
+        #bd = H2BladeData(r"C:\mmpe\HAWC2\models\NREL 5MW reference wind turbine_play/htc/NREL_5MW_reference_wind_turbine_heli_step.htc", "../")
+        if 1:
+            bd.plot_geometry()
+ 
+        
+
+
+    
+
+if __name__ == "__main__":
+    #import sys;sys.argv = ['', 'Test.testBladedata']
+    unittest.main()