diff --git a/wetb/hawc2/bladeinfo.py b/wetb/hawc2/bladeinfo.py
index 12c1d4aa556aa99c6facaabe0a2609d1711a98ce..e89bd9ab10c5d3619f38b442870bf0e115250cd8 100644
--- a/wetb/hawc2/bladeinfo.py
+++ b/wetb/hawc2/bladeinfo.py
@@ -14,53 +14,54 @@ from wetb.hawc2.st_file import StFile
-class BladeInfo(object):
- def twist(self, radius=None):
- """Aerodynamic twist [deg]. Negative when leading edge is twisted towards wind(opposite to normal definition)\n
+# class BladeInfo(object):
+# def twist(self, radius=None):
+# """Aerodynamic twist [deg]. Negative when leading edge is twisted towards wind(opposite to normal definition)\n
+#
+# Parameters
+# ---------
+# radius : float or array_like, optional
+# radius [m] of interest
+# If None (default) the twist of all points are returned
+#
+# Returns
+# -------
+# twist [deg] : float
+# """
+# raise NotImplementedError()
+#
+# def chord(self, radius=None):
+# """Aerodynamic chord
+#
+# Parameters
+# ---------
+# radius : float or array_like, optional
+# radius [m] of interest
+# If None (default) the twist of all points are returned
+#
+# Returns
+# -------
+# chord [m] : float
+# """
+# raise NotImplementedError()
+#
+#
+#
+# def c2nd(self, radius_nd):
+# """Return center line position
+#
+# Parameters
+# ---------
+# radius_nd : float or array_like, optional
+# non dimensional radius
+#
+# Returns
+# -------
+# x,y,z : float
+# """
- Parameters
- ---------
- radius : float or array_like, optional
- radius [m] of interest
- If None (default) the twist of all points are returned
-
- Returns
- -------
- twist [deg] : float
- """
- raise NotImplementedError()
-
- def chord(self, radius=None):
- """Aerodynamic chord
-
- Parameters
- ---------
- radius : float or array_like, optional
- radius [m] of interest
- If None (default) the twist of all points are returned
-
- Returns
- -------
- chord [m] : float
- """
- raise NotImplementedError()
-
-
- def c2nd(self, radius_nd):
- """Return center line position
-
- Parameters
- ---------
- radius_nd : float or array_like, optional
- non dimensional radius
-
- Returns
- -------
- x,y,z : float
- """
-
-class H2BladeInfo(BladeInfo, PCFile, AtTimeFile):
+class H2aeroBladeInfo(PCFile, AtTimeFile):
"""Provide HAWC2 info about a blade
From AE file:
@@ -85,45 +86,64 @@ class H2BladeInfo(BladeInfo, PCFile, AtTimeFile):
"""
- def __init__(self, htcfile, ae_filename=None, pc_filename=None, at_time_filename=None, st_filename=None, blade_name=None):
+ def __init__(self, htcfile, ae_filename=None, pc_filename=None, at_time_filename=None, blade_name=None):
if isinstance(htcfile, str):
assert htcfile.lower().endswith('.htc')
htcfile = HTCFile(htcfile)
-
+ self.htcfile = htcfile
blade_name = blade_name or htcfile.aero.link[2]
- s = htcfile.new_htc_structure
+ #s = htcfile.new_htc_structure
at_time_filename = at_time_filename or ("output_at_time" in htcfile and os.path.join(htcfile.modelpath, htcfile.output_at_time.filename[0] + ".dat"))
pc_filename = pc_filename or os.path.join(htcfile.modelpath, htcfile.aero.pc_filename[0])
ae_filename = ae_filename or os.path.join(htcfile.modelpath, htcfile.aero.ae_filename[0])
- mainbodies = [s[k] for k in s.keys() if s[k].name_ == "main_body"]
- mainbody_blade = [mb for mb in mainbodies if mb.name[0] == blade_name][0]
- st_filename = st_filename or os.path.join(htcfile.modelpath, mainbody_blade.timoschenko_input.filename[0])
+ #mainbodies = [s[k] for k in s.keys() if s[k].name_ == "main_body"]
+ #self.mainbody_blade = [mb for mb in mainbodies if mb.name[0] == blade_name][0]
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 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.c2def = np.array([v.values[1:5] for v in mainbody_blade.c2_def if v.name_ == "sec"])
+ self.curved_length = self.radius_curved_ac()[-1]
+ else:
+ raise NotImplementedError
+ #z_nd = (np.cos(np.linspace(np.pi, np.pi*2,len(curved_length)-1))+1)/2
+ #self.curved_length = np.cumsum(np.sqrt(np.sum(np.diff(self.c2def[:, :3], 1, 0) ** 2, 1)))[-1]
- #self.c2nd = lambda x : interp1d(self.c2def[:, 2] / self.c2def[-1, 2], self.c2def[:, :3], axis=0, kind='cubic')(np.max([np.min([x, np.ones_like(x)], 0), np.zeros_like(x) + self.c2def[0, 2] / self.c2def[-1, 2]], 0))
- x, y, z, twist = self.hawc2_splines()
- def interpolate(r):
- r = (np.max([np.min([r, np.ones_like(r)], 0), np.zeros_like(r) + self.c2def[0, 2] / self.c2def[-1, 2]], 0))
- return np.array([np.interp(r, np.array(z) / z[-1], xyz) for xyz in [x,y,z, twist]]).T
-
- self.c2nd = interpolate #lambda r : interp1d(np.array(z) / z[-1], np.array([x, y, z]).T, axis=0, kind=1)(np.max([np.min([r, np.ones_like(r)], 0), np.zeros_like(r) + self.c2def[0, 2] / self.c2def[-1, 2]], 0))
+
+ self.hawc2_splines_data = self.hawc2_splines()
+
+
+ @property
+ def c2def(self):
+ if not hasattr(self, "_c2def"):
+ self._c2def = np.array([v.values[1:5] for v in self.htcfile.blade_c2_def if v.name_ == "sec"])
+ return self._c2def
+
+
+ def c2nd(self, l_nd, curved_length=True):
+ curve_l_nd, x, y, z, twist = self.hawc2_splines_data
+ if curved_length:
+ l_nd = (np.max([np.min([l_nd, np.ones_like(l_nd)], 0), np.zeros_like(l_nd) + self.c2def[0, 2] / self.c2def[-1, 2]], 0))
+ return np.array([np.interp(l_nd, curve_l_nd, xyz) for xyz in [x,y,z, twist]]).T
+ else:
+ l_nd = (np.max([np.min([l_nd, np.ones_like(l_nd)], 0), np.zeros_like(l_nd)], 0))
+ return np.array([np.interp(l_nd, z/z[-1], xyz) for xyz in [x,y,z, twist]]).T
+
+ def c2(self, l, curved_length=True):
+ if curved_length:
+ L = self.curved_length
+ else:
+ L = self.c2def[-1,2]
+ return self.c2nd(l/L, curved_length)
+
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]
+ curve_l = np.r_[0, np.cumsum(np.sqrt(np.sum(np.diff(self.c2def[:, :3], 1, 0) ** 2, 1)))]
+ curve_l_nd = curve_l / curve_l[-1]
def akima(x, y):
n = len(x)
@@ -166,10 +186,11 @@ class H2BladeInfo(BladeInfo, PCFile, AtTimeFile):
z = np.linspace(0, s[i + 1] - s[i ], 10, endpoint=i >= co.shape[0] - 2)
x.extend(s[i] + z)
y.extend(p(z))
- return y
+ return x,y
- x, y, z, twist = [coef2spline(curve_z_nd, akima(curve_z_nd, self.c2def[:, i])) for i in range(4)]
- return x, y, z, twist
+ x, y, z, twist = [coef2spline(curve_l_nd, akima(curve_l_nd, self.c2def[:, i]))[1] for i in range(4)]
+ curve_l_nd = coef2spline(curve_l_nd, akima(curve_l_nd, self.c2def[:, 0]))[0]
+ return curve_l_nd, x, y, z, twist
def xyztwist(self, l=None, curved_length=False):
"""Return splined x,y,z and twist
@@ -188,44 +209,120 @@ class H2BladeInfo(BladeInfo, PCFile, AtTimeFile):
x,y,z,twist
"""
if l is None:
- return self.c2def[:, 3]
+ return self.c2def[:, :3]
else:
- r_nd = np.linspace(0,1,100)
+ r_nd = np.linspace(0,1,1000)
if curved_length:
- curved_length = np.cumsum(np.sqrt((np.diff(self.c2nd(np.linspace(0,1,100)),1,0)[:,:3]**2).sum(1)))
- assert np.all(l>=curved_length[0]) and np.all(l<=curved_length[-1])
- return self.c2nd(r_nd[np.argmin(np.abs(curved_length-l))+1])
+ #curved_length = np.cumsum(np.sqrt((np.diff(self.c2nd(np.linspace(0,1,100)),1,0)[:,:3]**2).sum(1)))
+ return self.c2nd(l/self.radius_curved_ac()[-1])
+# z_nd = (np.cos(np.linspace(np.pi, np.pi*2,len(curved_length)-1))+1)/2
+# assert np.all(l>=curved_length[0]) and np.all(l<=curved_length[-1])
+# return self.c2nd(r_nd[np.argmin(np.abs(curved_length-l))+1])
else:
assert np.all(l>=self.c2def[0,2]) and np.all(l<=self.c2def[-1,2])
return self.c2nd(l/self.c2def[-1, 2])
-
-
-class H2aeroBladeInfo(H2BladeInfo):
-
- def __init__(self, at_time_filename, ae_filename, pc_filename, htc_filename):
- """
- at_time_filename: file name of at time file containing twist and chord data
- """
- PCFile.__init__(self, pc_filename, ae_filename)
- blade_radius = self.ae_sets[1][-1,0]
- AtTimeFile.__init__(self, at_time_filename, blade_radius)
-
- assert('twist' in self.attribute_names)
- htcfile = HTCFile(htc_filename)
-
-
- self.c2def = np.array([v.values[1:5] for v in htcfile.blade_c2_def if v.name_ == "sec"])
- #self._c2nd = interp1d(self.c2def[:, 2] / self.c2def[-1, 2], self.c2def[:, :3], axis=0, kind='cubic')
-
- ac_radii = self.ac_radius()
- c2_axis_length = np.r_[0, np.cumsum(np.sqrt((np.diff(self.c2def[:, :3], 1, 0) ** 2).sum(1)))]
- self._c2nd = interp1d(c2_axis_length / c2_axis_length[-1], self.c2def[:, :3], axis=0, kind='cubic')
- #self._c2nd = interp1d(self.c2def[:,2]/self.c2def[-1,2], self.c2def[:, :3], axis=0, kind='cubic')
-
- def c2nd(self, r_nd):
- r_nd_min = np.zeros_like(r_nd) + self.c2def[0, 2] / self.c2def[-1, 2]
- r_nd_max = np.ones_like(r_nd)
- r_nd = np.max([np.min([r_nd, r_nd_max], 0), r_nd_min], 0)
- return self._c2nd(r_nd)
+
+
+class H2BladeInfo(H2aeroBladeInfo):
+ """Provide HAWC2 info about a blade
+
+ From AE file:
+ - chord(radius=None, set_nr=1):
+ - thickness(radius=None, set_nr=1)
+ - radius_ae(radius=None, set_nr=1)
+
+ From PC file
+ - CL(radius, alpha, ae_set_nr=1)
+ - CD(radius, alpha, ae_set_nr=1)
+ - CM(radius, alpha, ae_set_nr=1)
+ From at_time_filename
+ - attribute_names
+ - xxx(radius=None, curved_length=None) # xxx for each attribute name
+ - radius_curved_ac(radius=None) # Curved length of nearest/all aerodynamic calculation points
+
+ From ST file
+ - radius_st(radius=None, mset=1, set=1)
+ - xxx(radius=None, mset=1, set=1) # xxx for each of r, m, x_cg,y_cg, ri_x, ri_y, xs, ys, E, G, Ix, Iy, K, kx, ky, A, pitch, xe, ye
+ with template
+
+ """
+ def __init__(self, htcfile, ae_filename=None, pc_filename=None, at_time_filename=None, st_filename=None, blade_name=None):
+ if isinstance(htcfile, str):
+ htcfile = HTCFile(htcfile)
+ s = htcfile.new_htc_structure
+# at_time_filename = at_time_filename or ("output_at_time" in htcfile and os.path.join(htcfile.modelpath, htcfile.output_at_time.filename[0] + ".dat"))
+# pc_filename = pc_filename or os.path.join(htcfile.modelpath, htcfile.aero.pc_filename[0])
+# ae_filename = ae_filename or os.path.join(htcfile.modelpath, htcfile.aero.ae_filename[0])
+#
+ mainbodies = [s[k] for k in s.keys() if s[k].name_ == "main_body"]
+ if blade_name is None:
+ blade_name = htcfile.aero.link[2]
+ self.mainbody_blade = [mb for mb in mainbodies if mb.name[0] == blade_name][0]
+ H2aeroBladeInfo.__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):
+#
+# if isinstance(htcfile, str):
+# assert htcfile.lower().endswith('.htc')
+# htcfile = HTCFile(htcfile)
+#
+# blade_name = blade_name or htcfile.aero.link[2]
+ st_filename = st_filename or os.path.join(htcfile.modelpath, self.mainbody_blade.timoschenko_input.filename[0])
+#
+# 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 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]
+# else:
+# raise NotImplementedError
+# #z_nd = (np.cos(np.linspace(np.pi, np.pi*2,len(curved_length)-1))+1)/2
+# #self.curved_length = np.cumsum(np.sqrt(np.sum(np.diff(self.c2def[:, :3], 1, 0) ** 2, 1)))[-1]
+#
+# 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"):
+ self._c2def = np.array([v.values[1:5] for v in self.mainbody_blade.c2_def if v.name_ == "sec"])
+ return self._c2def
+
+#
+# class H2aeroBladeInfo(H2BladeInfo):
+#
+# def __init__(self, at_time_filename, ae_filename, pc_filename, htc_filename):
+# """
+# at_time_filename: file name of at time file containing twist and chord data
+# """
+# PCFile.__init__(self, pc_filename, ae_filename)
+# blade_radius = self.ae_sets[1][-1,0]
+# AtTimeFile.__init__(self, at_time_filename, blade_radius)
+#
+# assert('twist' in self.attribute_names)
+# htcfile = HTCFile(htc_filename)
+#
+#
+# self.c2def = np.array([v.values[1:5] for v in htcfile.blade_c2_def if v.name_ == "sec"])
+# #self._c2nd = interp1d(self.c2def[:, 2] / self.c2def[-1, 2], self.c2def[:, :3], axis=0, kind='cubic')
+#
+# ac_radii = self.radius_curved_ac()
+# self.hawc2_splines_data = self.hawc2_splines()
+# self.curved_length = np.cumsum(np.sqrt(np.sum(np.diff(self.c2def[:, :3], 1, 0) ** 2, 1)))[-1]
+#
+# # c2_axis_length = np.r_[0, np.cumsum(np.sqrt((np.diff(self.c2def[:, :3], 1, 0) ** 2).sum(1)))]
+# # self._c2nd = interp1d(c2_axis_length / c2_axis_length[-1], self.c2def[:, :3], axis=0, kind='cubic')
+# #self._c2nd = interp1d(self.c2def[:,2]/self.c2def[-1,2], self.c2def[:, :3], axis=0, kind='cubic')
+# #
+# # def c2nd(self, r_nd):
+# # r_nd_min = np.zeros_like(r_nd) + self.c2def[0, 2] / self.c2def[-1, 2]
+# # r_nd_max = np.ones_like(r_nd)
+# # r_nd = np.max([np.min([r_nd, r_nd_max], 0), r_nd_min], 0)
+# # return self._c2nd(r_nd)
+#