diff --git a/docs/source/bladed/Bladed_airfoil.png b/docs/source/bladed/Bladed_airfoil.png
new file mode 100644
index 0000000000000000000000000000000000000000..25bda289fe7e84f4518224469f4d3451a829c71d
Binary files /dev/null and b/docs/source/bladed/Bladed_airfoil.png differ
diff --git a/docs/source/bladed/Bladed_st_centers.png b/docs/source/bladed/Bladed_st_centers.png
new file mode 100644
index 0000000000000000000000000000000000000000..30772804ba8f5b09a8e989b4f9b04077af342aa1
Binary files /dev/null and b/docs/source/bladed/Bladed_st_centers.png differ
diff --git a/docs/source/bladed/Bladed_turbine_coord.png b/docs/source/bladed/Bladed_turbine_coord.png
new file mode 100644
index 0000000000000000000000000000000000000000..e9a897791ab413d694824edc4f08885338ad9259
Binary files /dev/null and b/docs/source/bladed/Bladed_turbine_coord.png differ
diff --git a/docs/source/bladed/HAWC2_c2_def_ccord.png b/docs/source/bladed/HAWC2_c2_def_ccord.png
new file mode 100644
index 0000000000000000000000000000000000000000..f940380b8f03988a8a653763efd800f9c7ca0722
Binary files /dev/null and b/docs/source/bladed/HAWC2_c2_def_ccord.png differ
diff --git a/docs/source/bladed/HAWC2_st_centers.png b/docs/source/bladed/HAWC2_st_centers.png
new file mode 100644
index 0000000000000000000000000000000000000000..1d0430f2be5ac61be4f3b8308712411049344d67
Binary files /dev/null and b/docs/source/bladed/HAWC2_st_centers.png differ
diff --git a/docs/source/bladed/HAWC2_turbine_coord.png b/docs/source/bladed/HAWC2_turbine_coord.png
new file mode 100644
index 0000000000000000000000000000000000000000..34e45d8e99dbb916aa3c564c721427e8cd4fabaf
Binary files /dev/null and b/docs/source/bladed/HAWC2_turbine_coord.png differ
diff --git a/docs/source/bladed/bladed2hawc.rst b/docs/source/bladed/bladed2hawc.rst
new file mode 100644
index 0000000000000000000000000000000000000000..20f4faaa4d1e51bd10e94010180a6afcd3f4e325
--- /dev/null
+++ b/docs/source/bladed/bladed2hawc.rst
@@ -0,0 +1,185 @@
+********************************
+BLADED to HAWC2 model conversion
+********************************
+
+**BLADED** is a tool developed by DNV, and ``wetb`` contains a basic reader that
+can interpret a BLADED project file (.prj) and convert the data structures into
+an equivalent set of **HAWC2** input files (htc, st, ae, pc).
+
+
+================================
+API and workflow
+================================
+
+---------------------------------
+Read BLADED project file
+---------------------------------
+
+The input from a Bladed model is given in a project file (*prj*) and is formatted
+in the XML definition. Within this file the core BLADED inputs are given in the
+CDATA field and which according to the XML standards refer to Character Data,
+see also `here <https://www.w3resource.com/xml/CDATA-sections.php>`_ and `here
+<https://en.wikipedia.org/wiki/CDATA>`_ for a more detailed technical description.
+Within the CDATA field a custom BLADED type markup definition is used. The ``wetb``
+reader includes an interpreter for this custom BLADED markup and outputs this
+data structure to the user as a nested dictionary, and which allows easy access
+to all the different key/value pairs defined in the BLADED CDATA field. This
+nested dictionary structure is accessible via the variable
+``wetb.bladed.readprj.ReadBladedProject.bd``.
+
+Outside the CDATA field the XML data is parsed via the ``lxml`` Python package,
+and the entire XML object ``lxml.objectify.fromstring(str)`` is exposed to the user
+via the ``wetb.bladed.readprj.ReadBladedProject.xmlroot`` variable.
+
+
+In ``wetb` the
+project file can be read using the following::
+
+ from wetb.bladed.readprj import ReadBladedProject
+ prj = ReadBladedProject('my_bladed_project_file.prj')
+ # XML object tree (and lxml data object)
+ prj.xmlroot
+ # CDATA field as nested dictionary
+ prj.bd
+
+
+A convenience function ``wetb.bladed.readprj.ReadBladedProject.get_key`` is also
+available to extract a ``numpy`` array, for example as follows::
+
+ from wetb.bladed.readprj import ReadBladedProject
+ prj = ReadBladedProject('my_bladed_project_file.prj')
+ # the following keys contain a data array
+ data_arr = prj.get_key('BSTIFFMB', 'EIFLAP')
+ # which is the similar calling (but excludes data conversion to int/float)
+ prj.bd['BSTIFFMB']['EIFLAP']
+ # if a key contains a section with sub-keys, a dictionary is returned:
+ data_dict = prj.get_key('BSTIFFMB')
+ # and is similar to
+ prj.bd['BSTIFFMB']
+
+
+---------------------------------
+Convert BLADED project file to HAWC2
+---------------------------------
+
+The class ``wetb.bladed.prj2hawc.Convert2Hawc`` will convert a Bladed project
+file into a set of HAWC2 input files. This process assumes that a standard
+3 bladed upwind turbine configuration is used, and a generic HAWC2 `htc template
+file <https://gitlab.windenergy.dtu.dk/toolbox/WindEnergyToolbox/-/blob/master/wetb/bladed/template.htc>`_
+serves as the starting point for the conversion process. Note that the converter
+class ``wetb.bladed.prj2hawc.Convert2Hawc`` inherits from
+``wetb.bladed.readprj.ReadBladedProject``, and hence will first read the Bladed
+project file::
+
+ import os
+ import urllib.request
+ import shutil
+ from wetb.bladed.prj2hawc import Convert2Hawc
+ # download the htc template file,
+ url = 'https://gitlab.windenergy.dtu.dk/toolbox/WindEnergyToolbox/-/'
+ url += 'raw/master/wetb/bladed/template.htc'
+ fname_tmpl = '/path/to/my/template.htc'
+ with urllib.request.urlopen(url) as response, open(fname_tmpl, 'wb') as fout:
+ shutil.copyfileobj(response, fout)
+ # file name of the Bladed project file
+ fprj = '/path/to/my/bladed.prj'
+ # initiate the project converter object, will first read the prj file (see above)
+ prj = Convert2Hawc(fprj)
+ # convert to the HAWC2 formats: htc, ae, pc and st files
+ fname_htc = fname_tmpl.replace('template.htc', 'bladed-to-hawc2.htc')
+ prj.convert(fname_htc)
+ # note that the ae, pc and st files will be saved in the same directory as
+ # fname_htc, with .htc replaced by .ae, .pc, .st etc.
+
+
+================================
+Theoretical background
+================================
+
+---------------------------------
+Coordinate systems
+---------------------------------
+
+The coordinate systems used in both codes are different. Note that in HAWC2
+coordinate systems of general bodies are defined by the users. However,
+the HAWC2 global and blade cross-section coordinate systems are predetermined,
+and they differ from Bladed in the following manner:
+
+* Global Bladed *X* axis is positive pointing downwind whereas in HAWC2 the
+ global *Y* direction is positive in the downwind direction.
+
+* The HAWC2 global *Z* direction is in the same direction as the gravity vector
+ (positive down), while it is the opposite in Bladed (positive up).
+
+* The Bladed cross-section coordinate system is rotated 90 degrees around *Z* wrt HAWC2.
+
+The figures below illustrate clearly the HAWC2 and Bladed coordinate systems.
+
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+Bladed coordinate systems
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+.. figure:: Bladed_turbine_coord.png
+ :width: 200
+ :alt: bladed_turbine_coord
+
+ Bladed coordinate system, rotor rotation and radius definition.
+
+.. figure:: Bladed_st_centers.png
+ :width: 400
+ :alt: bladed_st_coord
+
+ Bladed cross-section structural centers, half chord location and structural pitch definition.
+
+.. figure:: Bladed_airfoil.png
+ :width: 400
+ :alt: bladed_airfoil_coord
+
+ Bladed airfoil geometric positioning along the blade.
+
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+HAWC2 coordinate systems
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+.. figure:: HAWC2_turbine_coord.png
+ :width: 200
+ :alt: hawc2_turbine_coord
+
+ HAWC2 global, tower, shaft, hub, blade-i and meteorological coordinate systems.
+
+.. figure:: HAWC2_st_centers.png
+ :width: 400
+ :alt: hawc2_st_ccord
+
+ HAWC2 cross-section structural centers, half chord location and structural pitch definition.
+
+.. figure:: HAWC2_c2_def_ccord.png
+ :width: 400
+ :alt: hawc2_airfoil_coord
+
+ HAWC2 airfoil positioning in blade body coordinates and aerodynamic pitch given as in the *htc* file *c2_def* section.
+
+
+---------------------------------
+Cross-sectional parameters
+---------------------------------
+
+Bladed uses isotropic material definitions for all bodies and for the HAWC2
+conversion the same isotropic assumption is used. Since the HAWC2 *st* file definition
+splits the Young's (*E*) and shear modulus (*G*) from the actual stiffness terms, and Bladed
+defines the actual stiffness values (meaning the product of *EI* etc), the
+corresponding HAWC2 *st* input simply assumes a value for *E* and *G*, and
+specifies the inertia such that the product (i.e. stiffness) is correct.
+
+Bladed defines a mass polar moment of inertia in combination with the ratio
+between the mass radii of gyration around the airfoil's center of mass, while
+in HAWC2 the radii of gyration in *X* and *Y* direction are given wrt the elastic
+center.
+
+
+---------------------------------
+References
+---------------------------------
+
+[1] `HAWC2 User Manual v12.8 <http://tools.windenergy.dtu.dk/HAWC2/manual/>`_
+
+[2] Bladed 4.6 Manual
+
diff --git a/docs/source/index.rst b/docs/source/index.rst
index d13f3e8dd93bb6bc333145efbce889d7dec93ec0..b3c1c44073bb7bd66edc49d0bb4328d039931319 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -13,7 +13,7 @@ tab on the left to install the code and get started.
Source code repository (and issue tracker):
https://gitlab.windenergy.dtu.dk/toolbox/WindEnergyToolbox
-
+
License:
GPLv3_
@@ -25,4 +25,6 @@ Contents:
installation
fatigue_tools/fatigue
- hawc2/hawc2
\ No newline at end of file
+ hawc2/hawc2
+ bladed/bladed2hawc
+
diff --git a/wetb/bladed/prj2hawc.py b/wetb/bladed/prj2hawc.py
new file mode 100644
index 0000000000000000000000000000000000000000..803930e730af7afd48f41c3e02083807120a0c1c
--- /dev/null
+++ b/wetb/bladed/prj2hawc.py
@@ -0,0 +1,847 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Jun 23 11:20:28 2021
+
+@author: dave
+"""
+
+import os
+from os.path import join as pjoin
+
+import numpy as np
+import pandas as pd
+
+# from wetb.prepost import misc
+from wetb.hawc2 import (HTCFile, AEFile, PCFile, StFile)
+
+from readprj import ReadBladedProject
+
+
+# TODO: move to misc?
+def get_circular_inertia(d1, t, m):
+ d2 = d1 - 2*t
+ A = np.pi/4*(d1**2-d2**2)
+ I_g = np.pi/64*(d1**4-d2**4)
+ I_m = I_g/A*m
+ return A,I_g,I_m
+
+
+# TODO: move to HTCFile?
+def add_c2_arr_htc(htc, c2_arr, body):
+ """add a c2_def in the form of an array to a htc c2_def section. Needs to
+ be an existing body, and the current c2_def will be overwritten.
+ """
+ nr_nodes = c2_arr.shape[0]
+ c2_def = htc['new_htc_structure'].get_subsection_by_name(body)['c2_def']
+
+ keys = ['sec'] + [f'sec__{i}' for i in range(2, nr_nodes+1)]
+
+ c2_def['nsec'].values[0] = nr_nodes
+ for i, key in enumerate(keys):
+ # control the formatting, otherwise it looks very messy in the htc file
+ # this we do by converting the values to strings
+ vals_fmt = [f'{k: 15.07f}' for k in c2_arr[i,:].tolist()]
+ # node nr is stripped from leading spaces, so it doesn't help adding them
+ values = [i+1] + vals_fmt
+ if key in c2_def:
+ c2_def[key].values[:] = values[:]
+ else:
+ c2_def.add_line(name=key, values=values)
+
+ return htc
+
+
+# TODO: move to HTCFile?
+def get_c2def(htc, body_name):
+ """
+ """
+
+ # select the c2def section from the shaft
+ c2_def = htc['new_htc_structure'].get_subsection_by_name(body_name)['c2_def']
+ # safe all c2_def positions in one array
+ c2_arr = np.zeros((len(c2_def.keys())-1,4))
+ # iterate over all nsec lines, but ignore the line that gives the number of points
+ for i, sec in enumerate(c2_def.keys()[1:]):
+ c2_arr[i,:] = c2_def[sec].values[1:]
+
+ return pd.DataFrame(c2_arr, columns=['x', 'y', 'z', 'twist'])
+
+
+class Convert2Hawc(ReadBladedProject):
+ """Based on the BLADED user manual v4.8
+ """
+
+ def __init__(self, fname):
+ super().__init__(fname)
+
+ def get_pc(self, setname_filter=None):
+ # ---------------------------------------------------------------------
+ # we ignore the airfoils blocks at the start of the ADAT sub-section
+ # ---------------------------------------------------------------------
+ # nr_airfoils = prj.get_key('ADAT', 'NFOILS')[0,0]
+ # airfoil_keys = ['NFOIL', 'FTYPE', 'NTHICK', 'NREYN']
+ # airfoils = {}
+ # for k in range(nr_airfoils):
+ # key_app = ''
+ # if k > 0:
+ # key_app = f'__{k+1}'
+ # airfoils[f'NFOIL_{k+1}'] = {}
+ # for key in airfoil_keys:
+ # ---------------------------------------------------------------------
+
+ # total number of pc sets
+ nsets = self.get_key('ADAT', 'NSETS')[0,0]
+ key_app = [''] + [f'__{k}' for k in range(2, nsets+1)]
+
+ # FIXME: are all sets always named consistantly? Meaning only the
+ # thickness changes in the setname, and the rest remains the same?
+ # as a work-around: ignore sets that contains setname_filter in SETNAME
+ if setname_filter is not None:
+ drop = []
+ for i, ka in enumerate(key_app):
+ setname = self.get_key('ADAT', f'SETNAME{ka}')[0,0]
+ if setname.find(setname_filter) > -1:
+ drop.append(ka)
+ key_app = set(key_app) - set(drop)
+ # update nsets
+ nsets = len(key_app)
+
+ # initiate PCFile objects
+ pc = PCFile()
+ tc_arr = np.ndarray((nsets), dtype=np.float32)
+ pc_lst = []
+
+ for i, ka in enumerate(key_app):
+
+ # and assure cm is around quarter chord point
+ assert self.get_key('ADAT', 'XA')[0,0]==25
+
+ # prj.get_key('ADAT', f'SETNAME{ka}')
+ # prj.get_key('ADAT', f'SETNB{ka}')
+ # prj.get_key('ADAT', f'COMMENTS{ka}')
+ # prj.get_key('ADAT', f'DATE{ka}')
+ # prj.get_key('ADAT', f'REYN{ka}')
+ # prj.get_key('ADAT', f'DEPL{ka}')
+ # prj.get_key('ADAT', f'NALPHA{ka}')
+ # prj.get_key('ADAT', f'NVALS{ka}')
+
+ tc_arr[i] = self.get_key('ADAT', f'THICK{ka}')
+ shape = (self.get_key('ADAT', f'NALPHA{ka}')[0,0], 4)
+ pcs = np.ndarray(shape, dtype=np.float32)
+ pcs[:,0] = self.get_key('ADAT', f'ALPHA{ka}')[0,:]
+ pcs[:,1] = self.get_key('ADAT', f'CL{ka}')[0,:]
+ pcs[:,2] = self.get_key('ADAT', f'CD{ka}')[0,:]
+ pcs[:,3] = self.get_key('ADAT', f'CM{ka}')[0,:]
+ pc_lst.append(pcs)
+
+ # in BGEOM the airfoil number is referred to and not the thickness
+ # sort properly on thickness
+ isort = tc_arr.argsort()
+ pc.pc_sets = {1:(tc_arr[isort], [pc_lst[k] for k in isort])}
+
+ return pc
+
+ def get_ae(self):
+ # BGEOMMB defines for each element the start and end point. We only
+ # need half that of course.
+ chord = self.get_key('BGEOMMB', 'CHORD')[0,0::2]
+ # FIXME: is it radius or curved lenght
+ radius = self.get_key('BGEOMMB', 'RJ')[0,0::2]
+ thickness = self.get_key('BGEOMMB', 'BTHICK')[0,0::2]
+ # FIXME: multiple pc sets can be defined but are ignored for now.
+ pc_set_id = np.ones(chord.shape)
+
+ ae = AEFile()
+ ae.add_set(radius, chord, thickness, pc_set_id, set_id=None)
+
+ return ae
+
+ def get_blade_c2_def(self):
+
+ # BLADENAME
+ # DISTMODE RJ
+ # NBE 100
+ # FOIL
+ # MOVING
+
+ # CE_X
+ # CE_Y
+ diam = self.get_key('RCON', 'DIAM')[0,0]
+ cone = self.get_key('RCON', 'CONE')[0,0]
+ tilt = self.get_key('RCON', 'TILT')[0,0]
+
+ # in BLADED all is wrt LE, so we need chord
+ chord = self.get_key('BGEOMMB', 'CHORD')[0,0::2]
+
+ # c2_def = htc['new_htc_structure'].get_subsection_by_name('blade1')['c2_def']
+ twist = self.get_key('BGEOMMB', 'TWIST')[0,0::2]
+
+ # FIXME: what is the difference between RJ and DIST?
+ # radius can be given eiter as curved length or radius
+ # Distance: This can be entered as a distance along the blade or as a
+ # distance along the blade root Z-axis. Select the appropriate option
+ # at the base of the screen.
+ dist = self.get_key('BGEOMMB', 'DIST')[0,0::2]
+ rj = self.get_key('BGEOMMB', 'RJ')[0,0::2]
+ # reference frame is EC, that seems to be indicated on the plots of the
+ # blade rotor coordinates/pitch/axis tc (fig 4-1 and 4-2)
+ ref_x = self.get_key('BGEOMMB', 'REF_X')[0,0::2]
+ ref_y = self.get_key('BGEOMMB', 'REF_Y')[0,0::2]
+ # EC wrt to LE in percentage of chord
+ ce_x = self.get_key('BGEOMMB', 'CE_X')[0,0::2]
+ ce_y = self.get_key('BGEOMMB', 'CE_Y')[0,0::2]
+
+ # HAWC2 variables
+ theta = twist*-1
+ # theta_d = theta*180/np.pi
+ ea_x = (50-ce_y)*0.01*chord
+ ea_y = ce_x*0.01*chord
+
+ # safe all c2_def positions in one array
+ c2_arr = np.zeros((len(twist),4))
+ c2_arr[:,0] = -1*(ref_y + ((ea_x*np.cos(theta) - ea_y*np.sin(theta))))
+ c2_arr[:,1] = (ref_x - ((ea_x*np.sin(theta) + ea_y*np.cos(theta))))
+ # assume DIST is curved length
+ c2_arr[:,2] = rj #dist
+ c2_arr[:,3] = theta*180/np.pi
+
+ return c2_arr
+
+ def get_blade_st(self):
+ # BMASSMB
+ # BSTIFFMB
+
+ # in HAWC2 structural pitch is relative to the twist
+ twist = self.get_key('BGEOMMB', 'TWIST')[0,0::2]
+ dist = self.get_key('BGEOMMB', 'DIST')[0,0::2]
+ chord = self.get_key('BGEOMMB', 'CHORD')[0,0::2]
+ thickness = self.get_key('BGEOMMB', 'BTHICK')[0,0::2]
+
+ # CG location, in percentage chord wrt LE
+ cm_x = self.get_key('BMASSMB', 'CM_X')[0,0::2]
+ cm_y = self.get_key('BMASSMB', 'CM_Y')[0,0::2]
+ mass = self.get_key('BMASSMB', 'MASS')[0,0::2]
+ siner = self.get_key('BMASSMB', 'SINER')[0,0::2]
+
+ # Radii of gyration ratio: the radius of gyration of mass about y_m
+ # divided by the radius of gyration of mass about x_m. This defaults to
+ # the relative profile thickness but can be over-written by un-checking
+ # the “Use default radii of gyration ratio†checkbox.
+ rgratio = self.get_key('BMASSMB', 'RGRATIO')[0,0::2]
+ # FIXME: DEF_RGRATIO but which value belongs to which option?
+ # assume 0 means "default is off"
+ def_rgratio = self.get_key('BMASSMB', 'DEF_RGRATIO')[0,0]
+
+ # FIXME: BETA_M: figure 3-2 from v4.8: mass axis orientation, radians?
+ beta_m = self.get_key('BMASSMB', 'BETA_M')[0,0::2]
+ # Mass axis orientation: the orientation of the principle axis of inertia.
+ # This defaults to the orientation of aerodynamic twist, but can be
+ # over-written by un-checking the “Use default mass axis orientationâ€
+ # checkbox. (See diagram below)
+ # FIXME: DEF_BETA_M but which value belongs to which option?
+ # assume 0 means "default is off"
+ def_beta_m = self.get_key('BMASSMB', 'DEF_BETA_M')[0,0]
+
+ # DEF_RGRATIO 0
+ # DEF_BETA_M 0
+
+ # RADPOS
+ # PITCHPOS
+
+ # DEF_RGRATIO 0
+ # DEF_BETA_M 0
+ # ICEDBLADES 0, 0, 0
+ # ICEDENSITY 700
+ # TIPCHORD 0
+ # NPOINT 3
+ # RADPOS .2, 30.12345, 40.12345
+ # PITCHPOS .2, 30.1, 40.1
+ # PM_X 0, 0, 0
+ # PM_Y 20, 10, 20
+ # PMASS 1000, 30, 1
+
+ # elastic axis, in percentage chord wrt LE
+ ce_x = self.get_key('BGEOMMB', 'CE_X')[0,0::2]
+ ce_y = self.get_key('BGEOMMB', 'CE_Y')[0,0::2]
+
+ # 3.5 Blade stiffness distribution (page 19)
+ # The stiffness must be defined about the principal axis of inertia at
+ # each blade station (see 3.1). The stiffness is the product of Young’s
+ # Modulus for the material and the second moment of area for the xp or
+ # yp directions as appropriate. The principal axis orientation is defined
+ # as an input, and defaults to the aerodynamic twist. In this case it is
+ # assumed to be parallel and perpendicular to the chord line. If the
+ # principal axis orientation is different from the aerodynamic twist, click
+ # the Use default principal axis orientation to off. (see diagram below)
+ eiflap = self.get_key('BSTIFFMB', 'EIFLAP')[0,0::2]
+ eiedge = self.get_key('BSTIFFMB', 'EIEDGE')[0,0::2]
+ beta_s = self.get_key('BSTIFFMB', 'BETA_S')[0,0::2]
+ gj = self.get_key('BSTIFFMB', 'GJ')[0,0::2]
+ cs_x = self.get_key('BSTIFFMB', 'CS_X')[0,0::2]
+ cs_y = self.get_key('BSTIFFMB', 'CS_Y')[0,0::2]
+ gaflap = self.get_key('BSTIFFMB', 'GAFLAP')[0,0::2]
+ gaedge = self.get_key('BSTIFFMB', 'GAEDGE')[0,0::2]
+ # FIXME: DEF_BETA_S probably says what the frame of reference is
+ # assume 0 means "default is off"
+ # self.get_key('BMASSMB', 'DEF_BETA_S')[0,0::2]
+
+ # for plotting purposes
+ sel = beta_m < -1
+ beta_m[sel] += np.pi
+ sel = beta_s < -1
+ beta_s[sel] += np.pi
+
+ # calculate the curved length
+ # c2_arr = self.get_blade_c2_def()
+ # curved_len_arr = curved_length(c2_arr)
+
+ ea_x = (50-ce_y)*0.01*chord
+ ea_y = ce_x*0.01*chord
+ cg_x = (50-cm_y)*0.01*chord
+ cg_y = cm_x*0.01*chord
+
+ # mass moments
+ Ipm_y = siner + (mass*(cg_x-ea_x)**2)
+ ri_y = np.sqrt( (Ipm_y / (1+rgratio**2)) / mass )
+ ri_x = rgratio * ri_y
+ # in BLADED rgratio = y/x, but in HAWC2 it x=y so we swap them
+
+ # just choose E, G so the rest follows
+ E, G = 1e10, 1e09
+ # however, we don't know EA (probably because stiff?), so assume A is
+ # a full box based on chord length and maximum airfoil thickness
+ # this will make the extensional stiffness EA quite large, but the user
+ # would still end up with a physical meaningful value
+ A = chord*chord*thickness/100
+
+ starr = np.ndarray((len(dist),19))
+ starr[:,0] = dist
+ starr[:,1] = mass
+ starr[:,2] = cg_x
+ starr[:,3] = cg_y
+ # radius of gyration
+ starr[:,4] = ri_x
+ starr[:,5] = ri_y
+ # shear center
+ starr[:,6] = (50-cs_y)*0.01*chord
+ starr[:,7] = cs_x*0.01*chord
+ # E, G: choose
+ starr[:,8] = E
+ starr[:,9] = G
+ starr[:,10] = eiflap/E
+ starr[:,11] = eiedge/E
+ starr[:,12] = gj/G
+ # shear factor
+ starr[:,13] = gaflap/(G*A)
+ starr[:,14] = gaedge/(G*A)
+ starr[:,15] = A
+ starr[:,16] = -1*(beta_s - twist) # structural pitch
+ starr[:,17] = ea_x
+ starr[:,18] = ea_y
+
+ st = StFile()
+ st.main_data_sets = {1:{1:starr}}
+ # st.cols = wetb.hawc2.st_file.stc.split()
+ return st
+
+ def get_tower_c2_def_st(self):
+ """Ignores nodes that are not at the centerline
+
+ Returns
+ -------
+ c2_arr : TYPE
+ DESCRIPTION.
+ st : TYPE
+ DESCRIPTION.
+
+ """
+
+ t_n_e = int(self.get_key('TGEOM', 'NTE'))
+ t_n_n = int(self.get_key('TGEOM', 'NTS'))
+ t_el_mat = self.get_key('TGEOM', 'ELSTNS')
+ t_n_coord = self.get_key('TGEOM', 'TCLOCAL')
+ # tdiam = self.get_key('TGEOM', 'SEAL')
+ # tdiam = self.get_key('TGEOM', 'HYDRO')
+ towmgt = self.get_key('TGEOM', 'TOWMGT')
+ # MATERIAL, NOMATS
+ material = self.get_key('TMASS', 'MATERIAL')[0,0::2]
+ nomats = self.get_key('TMASS', 'NOMATS')
+
+ # Index of tower nodes (nodes at zero x-y positions)
+ t_a = np.where((abs(t_n_coord[:,0])+abs(t_n_coord[:,1])) == 0)[0]
+ t_nodes = t_n_coord[t_a]
+ # sort index according to z position of the nodes
+ n_ind = t_a[np.argsort(t_nodes[:,2])]
+ t_nodes = t_n_coord[n_ind]
+ # Index for tower elemetns
+ # Elemet indexes which does not contain both tower nodes
+ ind_not = [i for i,j in enumerate(t_el_mat) if not (j[0]-1) in n_ind
+ or not (j[1]-1) in n_ind]
+ # The elemetn indexes excluding the ind_not
+ ind_init = np.array([i for i in range(t_el_mat.shape[0]) if not i in ind_not])
+ sort_el = []
+ for i in n_ind[:-1]:
+ sort_el.append(np.where(t_el_mat[ind_init][:,0] == i+1)[0][0])
+ ind_el = ind_init[sort_el]
+
+ # safe all c2_def positions in one array
+ c2_arr = np.zeros((t_nodes.shape[0],4))
+ c2_arr[:,0:2] = t_nodes[:,:-1].copy()
+ c2_arr[:,2] = -t_nodes[:,-1] + t_nodes[0,-1] # It is minus for HAWC2 model
+ c2_arr[:,3] = 0.0
+ # element properties
+ t_d = self.get_key('TGEOM', 'TDIAM')[ind_el]
+ material = self.get_key('TMASS', 'MATERIAL')[0,0::2][ind_el]
+ t_thick = self.get_key('TMASS', 'WALLTHICK')[ind_el]
+ towmgt = self.get_key('TGEOM', 'TOWMGT')[ind_el]
+ # Mass
+ t_m = self.get_key('TMASS', 'TOWM')[ind_el]
+ t_p_m = self.get_key('TMASS', 'TOWPMI')[ind_el] # point mass for elements and nodes
+ t_flood = self.get_key('TMASS', 'FLOODED')[ind_el]
+ # material props from prj file
+ name_mat = np.unique(material)
+ mat_prop = {}
+ for i in name_mat:
+ mat_prop[i] = self.get_key('TMASS', "%s"%i) # rho, E, G
+
+ t_mat = np.array([mat_prop[i][0] for i in material])
+
+ if np.sum(t_thick[1:,0] - t_thick[:-1,1]) > 1e-9:
+ print('Diameters are not continous along the tower')
+ if (name_mat.shape[0] > 1):
+ print('There are %i materials in tower!' %name_mat.shape[0])
+ d = np.zeros(t_d.shape[0]+1)
+ t = np.zeros(t_thick.shape[0]+1)
+ mat = np.zeros((t_mat.shape[0]+1,t_mat.shape[1]))
+ d[1:] = t_d[:,1]
+ t[1:] = t_thick[:,1]
+ mat[1:,:] = t_mat.copy()
+ d[0] = t_d[0,1]
+ t[0] = t_thick[0,1]
+ mat[0,:] = t_mat[0,:].copy()
+
+ A,I_g,I_m = get_circular_inertia(d,t,mat[:,0])
+ # FIXME: MAYBE ANOTHER FIX FOR DIFFERENT MATERIALS
+ # return c2_arr
+ starr = np.zeros((c2_arr.shape[0],19))
+ starr[:,0] = -c2_arr[:,2] #
+ starr[:,1] = A*mat[:,0]
+ starr[:,2] = 0.0 # no cg offset
+ starr[:,3] = 0.0 # no cg offset
+ # radius of gyration
+ starr[:,4] = np.sqrt(I_m/A) #ri_x
+ starr[:,5] = np.sqrt(I_m/A) #ri_y
+ # shear center
+ starr[:,6] = 0.0 # no shear center offset
+ starr[:,7] = 0.0 # no shear center offset
+ # E, G: choose
+ starr[:,8] = mat[:,1]
+ starr[:,9] = mat[:,2]
+ starr[:,10] = I_g
+ starr[:,11] = I_g
+ starr[:,12] = I_g*2
+ # shear factor
+ starr[:,13] = 3/4 # for circular section check it again
+ starr[:,14] = 3/4 # for circular section check it again
+ starr[:,15] = A
+ starr[:,16] = 0.0 # structural pitch
+ starr[:,17] = 0.0
+ starr[:,18] = 0.0
+
+ st = StFile()
+ st.main_data_sets = {1:{1:starr}}
+
+ return c2_arr, st
+
+ def get_hub(self):
+ # Blade root length, in blade root coordinates (with tilt and cone)
+ # see figure 4-1 in Blade User Manual 4.8
+ root = self.get_key('RCON', 'ROOT')[0,0]
+
+ # Hub mass: the mass of the hub, including the spinner and any blade
+ # root section.
+ # Hub mass centre: the distance from the intersection of the shaft and
+ # blade axes to the centre of mass of the hub, in a direction measured
+ # away from the tower.
+ # Moments of inertia: the moment of inertia of the hub mass about the
+ # shaft axis must be defined.
+ # The inertia about an axis perpendicular to the shaft may also be
+ # entered with its origin about the hub centre of mass.
+
+ hubmas = self.get_key('RMASS', 'HUBMAS')[0,0]
+ hubine = self.get_key('RMASS', 'HUBINE')[0,0]
+ # FIXME: which axis is this? Would that translate to both xx and yy in H2?
+ hubine2 = self.get_key('RMASS', 'HUBINE2')[0,0]
+
+ # Enter the Spinner diameter. This is the diameter of any spinner or
+ # nose-cone, within which the blades themselves experience no
+ # aerodynamic forces.
+ spind = self.get_key('RCON', 'SPIND')[0,0]
+
+ cmass_hub = {'x':0, 'y':0, 'z':0, 'm':hubmas,
+ 'Ixx':0, 'Iyy':0, 'Izz':hubine}
+
+ return cmass_hub, root
+
+ def get_drivertain(self, len_shaft):
+
+ # TODO: brake
+ # MSTART BRAKE
+
+ # DTRAIN
+ # dtrain = self.get_key('DTRAIN')
+ ginert = self.get_key('DTRAIN', 'GINERT')[0,0] # generator inertia
+ gratio = self.get_key('DTRAIN', 'GRATIO')[0,0]
+ # The additional inertia of the high speed shaft may also be specified
+ # along with the inertia of the gearbox which is referred to the HSS.
+ gbxinert = self.get_key('DTRAIN', 'GBXINERT')[0,0] # gearbox inertia
+ # FIXME: brake position?
+ bpos = self.get_key('DTRAIN', 'BPOS')[0,0]
+ # LSS seems to have a DOF indicator in LSSDOF, and has 6 elements
+ # 4th element is torsion it seems
+ klss = self.get_key('DTRAIN', 'KLSS')[0,3]
+ # FIXME: what does this damping value mean? Simply C in the dyn system
+ # theta_dt_dt*I + C*theta_dt + K*theta = 0 ?
+ dlss = self.get_key('DTRAIN', 'DLSS')[0,3]
+ # HSS only has one DOF
+ khss = self.get_key('DTRAIN', 'KHSS')[0,0]
+ dhss = self.get_key('DTRAIN', 'DHSS')[0,0]
+
+ # Convert torsional stiffness definition into a beam stiffness
+ # that applies for the considered shaft length (c2_def)
+ G = 80e9 # typical steel value
+ Ip = klss*len_shaft/G
+ starr = np.zeros((2,19))
+ starr[:,0] = [0, len_shaft] # st file is always normalised
+ starr[:,1] = 0.5 # very light but not too light
+ starr[:,2] = 0.0 # no cg offset
+ starr[:,3] = 0.0 # no cg offset
+ # radius of gyration: very light but not too light
+ starr[:,4] = 0.001 #ri_x
+ starr[:,5] = 0.001 #ri_y
+ # shear center
+ starr[:,6] = 0.0 # no shear center offset
+ starr[:,7] = 0.0 # no shear center offset
+ # E, G: choose
+ starr[:,8] = 1e18
+ starr[:,9] = G
+ starr[:,10] = Ip/2 # area inertia's
+ starr[:,11] = Ip/2
+ starr[:,12] = Ip
+ # shear factor
+ starr[:,13] = 1000 # stiff in shear
+ starr[:,14] = 1000
+ starr[:,15] = 10 # cross-section area
+ starr[:,16] = 0.0 # structural pitch
+ starr[:,17] = 0.0 # elastic axis
+ starr[:,18] = 0.0
+ st = StFile()
+ st.main_data_sets = {1:{1:starr}}
+
+ # FIXME: INERTIATOCLUTCH??
+ inertiatoclutch = self.get_key('DTRAIN', 'INERTIATOCLUTCH')
+
+ # FIXME: what is the location of these inertia's?
+ # FIXME: what is the mass? It seems the mass is put elsewhere?
+
+ # cmass_gbx = {'Izz':Izz}
+ # cmass_gen = {'x':x, 'y':y, 'z':z, 'm':nacmas,
+ # 'Ixx':Ixx, 'Iyy':Iyy, 'Izz':Izz}
+
+ cmass_lss = {'x':0, 'y':0, 'z':0, 'm':0,'Ixx':0, 'Iyy':0, 'Izz':ginert}
+ # convert HSS inertia to LSS
+ cmass_hss = {'x':0, 'y':0, 'z':0, 'm':0,'Ixx':0, 'Iyy':0,
+ 'Izz':gbxinert*gratio*gratio}
+
+ return cmass_lss, cmass_hss, st
+
+ def get_nacelle(self):
+
+ # also get the nacell concentrated mass
+ nacmas = self.get_key('NMASS', 'NACMAS')[0,0]
+ # includes the nacelle structure and all the machinery within it.
+ # It does not include the rotor blades and hub. If the Direct Drive
+ # (see 4.1) option is selected, the mass of the generator (see 4.2)
+ # is also excluded.
+
+ # position of cg relative to tower top, nacelle mass coord system
+ # based on figure 4-7 of bladed manual 4.8
+ # FIXME: is this correct? z means SS, and at the edge of the nacelle?
+ nmx = self.get_key('NMASS', 'NMX')[0,0] # width?
+ nmy = self.get_key('NMASS', 'NMY')[0,0] # length?
+ nmz = self.get_key('NMASS', 'NMZ')[0,0] # height?
+ iyaw = self.get_key('NMASS', 'IYAW')[0,0]
+ inod = self.get_key('NMASS', 'INOD')[0,0]
+ iroll = self.get_key('NMASS', 'IROLL')[0,0]
+
+ # MSTART NGEOM
+ # NACW x.000
+ # NACL x.000
+ # NACH x.000
+ # NACZ 0
+ # NACCD x.0
+ # NAERO N
+
+ # height of the nacelle/towertop: the BLADED model does not have to
+ # have consistent shaft length with tilt angle arriving at the rotor
+ # center, so adjust towertop length in HAWC2 accordingly
+
+ # Tower height: from the ground or sea surface to the yaw bearing
+ # (only needed if the tower itself has not been defined)
+ towht = self.get_key('RCON', 'TOWHT')[0,0]
+
+ # Hub height: from the ground to the centre of the rotor
+ # (i.e. the intersection of the blade and shaft axes)
+ height = self.get_key('RCON', 'HEIGHT')[0,0]
+
+ # horizontal distance between rotor center and tower center line
+ ovrhng = self.get_key('RCON', 'OVRHNG')[0,0]
+ tilt = self.get_key('RCON', 'TILT')[0,0]
+
+ # FOR HAWC2, WE HAVE
+ len_shaft = ovrhng / np.cos(tilt)
+ len_towertop = height - towht - len_shaft*np.sin(tilt)
+ assert len_towertop > 0
+ # assert len_towertop + len_shaft*np.sin(tilt) + towht - height == 0
+
+ # concentrated mass in HAWC2
+ # 1 : nodenr
+ # 2-4: x, y, z offset
+ # 5 : mass
+ # 6-8: Ixx, Iyy, Izz
+ # 9-11: Ixy, Ixz, Iyz (optional)
+ # in HAWC2 we would be in the tower top coordinate system, and is
+ # usually the same as the global coordinate system
+ x = nmx
+ y = -nmy # move upwind wrt tower
+ z = -nmz # move upward
+ Ixx = inod
+ Iyy = iroll
+ Izz = iyaw
+
+ cmass_nacelle = {'x':x, 'y':y, 'z':z, 'm':nacmas,
+ 'Ixx':Ixx, 'Iyy':Iyy, 'Izz':Izz}
+
+ return cmass_nacelle, len_shaft, len_towertop
+
+ def get_towershadow_drag(self):
+ # TODO: tower shadow and tower aerodrag
+
+ # MSTART TOWSDW
+ # TSMODEL x
+ # TDCORR x
+ # MAXDEF x.x
+ # WIDTH x
+ # LREF x.x
+
+ return
+
+ def get_control(self):
+
+ # MSTART CONTROL
+
+ # MSTART PCOEFF
+ # TSRMIN x
+ # TSRMAX x
+ # TSRSTP x.x
+ # PITCH -x.0E-00
+ # PITCH_END -x.0E-00
+ # PITCH_STEP 0
+ # OMEGA x.654321
+
+ # MSTART IDLING
+
+ # AEROINFO
+
+ # pitch actuator rate limits
+ pitch_actuator = self.xmlroot.PitchActuator
+ if pitch_actuator.SetpointTrajectory.HasRateLimits:
+ # pitch_actuator.SetpointTrajectory.LowerRateLimit
+ # pitch_actuator.SetpointTrajectory.UpperRateLimit
+ print(pitch_actuator.SetpointTrajectory.UpperRateLimit*30/np.pi)
+
+ if pitch_actuator.SetpointTrajectory.HasAccelerationLimits:
+ # pitch_actuator.SetpointTrajectory.LowerAccelerationLimit
+ # pitch_actuator.SetpointTrajectory.UpperAccelerationLimit
+ print(pitch_actuator.SetpointTrajectory.UpperAccelerationLimit*30/np.pi)
+
+ # pitch_actuator.PositionResponse.values() # 1st or 2nd order
+ # pitch_actuator.PositionResponse.Frequency
+ # pitch_actuator.PositionResponse.Damping
+
+ # pitch_actuator.RateResponse.values() # 1st or 2nd order
+ # pitch_actuator.RateResponse.LagTime
+ # pitch_actuator.PositionResponse.Damping
+
+ # and more: bearing friction, etc
+
+ return
+
+ def convert(self, fname_htc):
+
+ fname_prj = os.path.basename(fname_htc).replace('.htc', '')
+
+ # for convenience, we start from an htc template that otherwise has the
+ # right structure, but has only 2 nodes for each body.
+ # assume it is in the same folder as where the HTC will be written to
+ basepath = os.path.dirname(fname_htc)
+ htc = HTCFile(pjoin(basepath, 'template.htc'))
+
+ # prj = self.PRJ_to_HAWC(pjoin(version, fname_prj))
+ # plot_bladed(prj)
+
+ # cone = prj.get_key('RCON', 'CONE')
+ # Lateral offset: the horizontal offset between the shaft and tower axes.
+
+ # Brake
+ # steps = prj.get_key('BRAKE', 'STEPS')
+ # TORQUE = prj.get_key('BRAKE', 'TORQUE')
+ # TIME = prj.get_key('BRAKE', 'TIME')
+
+ # -------------------------------------------------------------------------
+ # BLADES
+ # extract blade geometry, and add to c2_def section in the htc file
+ c2_blade = self.get_blade_c2_def()
+ htc = add_c2_arr_htc(htc, c2_blade, 'blade1')
+
+ # -------------------------------------------------------------------------
+ # NACELLE
+ cmass_nacelle, len_shaft, len_towertop = self.get_nacelle()
+ # set towertop length
+ # nacelle = htc['new_htc_structure'].get_subsection_by_name('towertop')
+ c2_tt = np.zeros((2,4))
+ c2_tt[1,2] = -len_towertop
+ htc = add_c2_arr_htc(htc, c2_tt, 'towertop')
+
+ # -------------------------------------------------------------------------
+ # attach nacelle mass/inertia to towertop, first node (yaw bearing)
+ towertop = htc['new_htc_structure'].get_subsection_by_name('towertop')
+ values = [v for k,v in cmass_nacelle.items()]
+ cmass = towertop.add_line('concentrated_mass', [1] + values)
+ key = cmass.location().split('/')[-1]
+ towertop[key].comments = 'nacelle mass (inc rotating) and inertia (non-rotating)'
+
+ # -------------------------------------------------------------------------
+ # set shaft length
+ c2_s = np.zeros((2,4))
+ c2_s[1,2] = len_shaft
+ htc = add_c2_arr_htc(htc, c2_s, 'shaft')
+
+ # -------------------------------------------------------------------------
+ # generator and gearbox inertia's, shaft torsional flexibility
+ cmass_lss, cmass_hss, st_shaft = self.get_drivertain(len_shaft)
+ st_shaft.save(pjoin(basepath, f'{fname_prj}_shaft.st'))
+
+ shaft = htc['new_htc_structure'].get_subsection_by_name('shaft')
+ # LSS inertia
+ values = [v for k,v in cmass_lss.items()]
+ cmass = shaft.add_line('concentrated_mass', [2] + values)
+ key = cmass.location().split('/')[-1]
+ shaft[key].comments = 'inertia LSS'
+ # Inertia of the HSS is already expressed wrt LSS
+ values = [v for k,v in cmass_hss.items()]
+ cmass = shaft.add_line('concentrated_mass', [2] + values)
+ gratio = self.get_key('DTRAIN', 'GRATIO')[0,0]
+ key = cmass.location().split('/')[-1]
+ shaft[key].comments = f'inertia HSS, expressed in LSS, GBR={gratio}'
+
+ # -------------------------------------------------------------------------
+ # set tilt angle
+ tilt = self.get_key('RCON', 'TILT')[0,0]
+ ori = htc['new_htc_structure'].orientation
+ # tilt is on the 2nd relative block, and set angle on 2nd set of eulerang
+ # check comments of the htc file:
+ # print(ori['relative__2']['mbdy2_eulerang__2'].comments)
+ ori['relative__2']['mbdy2_eulerang__2'].values = [tilt*180/np.pi, 0, 0]
+ ori['relative__2']['mbdy2_eulerang__2'].comments = 'tilt angle'
+
+ # coning angle
+ cone = self.get_key('RCON', 'CONE')[0,0]
+ # coning is on 3, 4 and 5 (all hub orientations)
+ for k in range(3,6):
+ # print(ori[f'relative__{k}']['mbdy2_eulerang__3'].comments)
+ ori[f'relative__{k}']['mbdy2_eulerang__3'].values = [-cone*180/np.pi, 0, 0]
+ ori[f'relative__{k}']['mbdy2_eulerang__3'].comments = 'cone angle'
+
+ # TODO: blade mounting sweep and cone as well??
+
+ # -------------------------------------------------------------------------
+ # HUB mass as concentrated mass on the shaft end
+ # FIXME: make sure we can get the hub root loads including these inertia's
+ cmass_hub, hublen = self.get_hub()
+ values = [v for k,v in cmass_hub.items()]
+ cmass = shaft.add_line('concentrated_mass', [2] + values)
+ key = cmass.location().split('/')[-1]
+ shaft[key].comments = 'Hub inertia and mass'
+
+ # -------------------------------------------------------------------------
+ # hub lenght
+ c2_hub = np.zeros((2,4))
+ c2_hub[1,2] = hublen
+ htc = add_c2_arr_htc(htc, c2_hub, 'hub1')
+
+ # -------------------------------------------------------------------------
+ # TOWER ST FILE and C2_DEF section
+ # WARNING: only includes nodes that are at centerline, side-elements are
+ # ignored
+ # TODO: compare between TGEOM/TMASS and TSTIFF approaches, it seems that
+ # the shear stiffness is not the same, see check_tower()
+ c2_tow, st = self.get_tower_c2_def_st()
+ st.save(pjoin(basepath, f'{fname_prj}_tower.st'))
+ htc = add_c2_arr_htc(htc, c2_tow, 'tower')
+
+ # hub height
+ htc.wind.center_pos0.values = [0, 0, -self.get_key('RCON', 'HEIGHT')[0,0]]
+
+ # -------------------------------------------------------------------------
+ # set all file names correct
+ tower = htc['new_htc_structure'].get_subsection_by_name('tower')
+ tower.timoschenko_input.filename.values[0] = f'data/{fname_prj}_tower.st'
+
+ shaft = htc['new_htc_structure'].get_subsection_by_name('shaft')
+ shaft.timoschenko_input.filename.values[0] = f'data/{fname_prj}_shaft.st'
+
+ blade1 = htc['new_htc_structure'].get_subsection_by_name('blade1')
+ blade1.timoschenko_input.filename.values[0] = f'data/{fname_prj}_blade.st'
+
+ for body_name in ['towertop', 'hub1']:
+ body = htc['new_htc_structure'].get_subsection_by_name(body_name)
+ body.timoschenko_input.filename.values[0] = 'data/template.st'
+
+ htc['aero'].ae_filename.values[0] = f'data/{fname_prj}.ae'
+ htc['aero'].pc_filename.values[0] = f'data/{fname_prj}.pc'
+
+ # -------------------------------------------------------------------------
+ # SAVE HTC FILE
+ htc.save(pjoin(basepath, f'{fname_htc}'))
+
+ # -------------------------------------------------------------------------
+ # other data files
+
+ # Blade st file
+ st = self.get_blade_st()
+ st.save(pjoin(basepath, f'{fname_prj}_blade.st'))
+
+ # extract profile coefficients and save to pc file
+ pc = self.get_pc(setname_filter='cleanVG')
+ pc.save(pjoin(basepath, f'{fname_prj}.pc'))
+
+ # extract aerodynamic layout, and safe to ae file
+ ae = self.get_ae()
+ ae.save(pjoin(basepath, f'{fname_prj}.ae'))
+
+ # # BLADE c2_def
+ # c2_def = htc['new_htc_structure'].get_subsection_by_name('blade1')['c2_def']
+ # # iterate over all nsec lines, but ignore the line that gives the number of points
+ # for i, sec in enumerate(c2_def.keys()[1:]):
+ # c2_arr[i,:] = c2_def[sec].values[1:]
+
+ # curvedlength = np.zeros(len(c2_arr[:,2]))
+ # for i in range(1,len(c2_arr[:,2])):
+ # curvedlength[i]=curvedlength[i-1]+np.linalg.norm(c2_arr[i,0:3]-c2_arr[i-1,0:3])
+
+ # c2_arr = np.zeros((len(c2_def.keys())-1,4))
+ # # iterate over all nsec lines, but ignore the line that gives the number of points
+ # for i, sec in enumerate(c2_def.keys()[1:]):
+ # c2_arr[i,:] = c2_def[sec].values[1:]
diff --git a/wetb/bladed/readprj.py b/wetb/bladed/readprj.py
new file mode 100644
index 0000000000000000000000000000000000000000..4d5291c627de31ca799c34a59fd880616a526285
--- /dev/null
+++ b/wetb/bladed/readprj.py
@@ -0,0 +1,344 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Jun 23 11:07:50 2021
+
+@author: dave
+"""
+
+# from os.path import join as pjoin
+
+import numpy as np
+from lxml import objectify#, etree)
+
+# from matplotlib import pyplot as plt
+# import pandas as pd
+
+# import wetb
+from wetb.prepost import misc
+# from wetb.hawc2 import (HTCFile, AEFile, PCFile, StFile)
+
+
+# TODO: this dictionary utility could go to a dict sub-class in misc?
+def find_next_unique_key(d, key):
+ k = 1
+ while True:
+ k += 1
+ if key + f'__{k}' not in d:
+ key = key + f'__{k}'
+ return key
+
+
+# TODO: this dictionary utility could go to a dict sub-class in misc?
+def add_new_unique_key(key, values, d):
+ """
+ Add key:values to dictionary d. If key already occurs, append __x
+ where x is the number of previous occurances of key(__x).
+
+ Parameters
+ ----------
+ key : str
+ DESCRIPTION.
+ values : list
+ DESCRIPTION.
+ d : dict
+ Dictionary.
+
+ Returns
+ -------
+ key : str
+ DESCRIPTION.
+ d : dict
+ DESCRIPTION.
+
+ """
+ if key in d:
+ key = find_next_unique_key(d, key)
+ d[key] = [values]
+ return key, d
+
+
+class ReadBladedProject:
+
+ def __init__(self, fname):
+
+ with open(fname, encoding='utf-8') as fobj:
+ xml_str = fobj.read().encode('utf-8')
+ self.bd, self.xmlroot = self.parse_bladeddata(xml_str)
+
+ self.set_attr_and_check()
+
+ # some things are just a little different
+ # TMASS has a list of materials and their properties
+ # get rid of the quotes
+ # tmp = [el.replace("'", '') for el in self.bd['TMASS']['MATERIAL'][0]]
+ # self.bd['TMASS']['MATERIAL'] = tmp
+ unique_mat = set(self.get_key('TMASS', 'MATERIAL').flatten().tolist())
+ self.tow_mat_prop = {k:self.get_key('TMASS', k) for k in unique_mat}
+ # material_props = {}
+ # for k in unique_mat:
+ # material_props[k.replace("'", '')] = self.get_key('TMASS', k)
+
+ def parse_bladeddata(self, xml_str):
+ """
+ The XML field BladedData contains what seems like the main core input
+ data for BLADED, and formatted in some structured way.
+
+ Parameters
+ ----------
+ xml_str : TYPE
+ DESCRIPTION.
+
+ Returns
+ -------
+ bd : TYPE
+ DESCRIPTION.
+ xmlroot : TYPE
+ DESCRIPTION.
+
+ """
+
+ # root = etree.fromstring(xml)
+ # elems = root.getchildren()
+ # bladeddata = elems[1].text
+
+ xmlroot = objectify.fromstring(xml_str)
+
+ # the non-xml formatted BLADED model is situated in a non-xml field
+ # called BladedData
+ bd = {}
+ mstart = None
+ for i, line in enumerate(xmlroot.BladedData.text.split('\n')):
+
+ # TODO: values embedded in double quotes (") can contain entire
+ # MSTART/MEND sub-sections (so embedded sections)
+
+ # split replace tabs with spaces, split on spaces, remove empty
+ linels = misc.remove_items(line.replace('\t', ' ').split(' '), '')
+ # commas can also be separators, in addition to spaces
+ linels2 = []
+ for k in linels:
+ linels2.extend(k.split(','))
+ linels = misc.remove_items(linels2, '')
+
+ # ignore empty lines
+ if len(linels) < 1:
+ continue
+
+ # entries can be numbers if the previous key contains multiple data points
+ try:
+ float(linels[0])
+ el0isnum = True
+ except ValueError:
+ el0isnum = False
+
+ # start of a sub-section that contains (non-unique) keys as well
+ if linels[0].upper().startswith('MSTART'):
+ mtag = linels[-1]
+ mstart = {}
+
+ # at the end of the sub-section add the sub-section to the main dict
+ elif linels[0].upper().startswith('MEND'):
+ # FIXME: implement MSTART sections embedded in double quoted values
+ try:
+ # if the section key is not unique, make it so by appending __X
+ if mtag in bd:
+ mtag = find_next_unique_key(bd, mtag)
+ bd[mtag] = mstart
+ except UnboundLocalError:
+ print('warning: ignored embedded mstart/mend section')
+ print(f'at line: {i+1}')
+ mstart = None
+
+ # if we are under a sub-section
+ elif mstart is not None:
+ # if the line contains a keyword
+ if not el0isnum:
+ tag, mstart = add_new_unique_key(linels[0], linels[1:], mstart)
+ # line is datapoint that needs to be added to key that occured before
+ else:
+ mstart[tag].append(linels)
+
+ # add numerical values to key that occured before
+ elif el0isnum:
+ bd[tag].append(linels)
+
+ else:
+ tag, bd = add_new_unique_key(linels[0], linels[1:], bd)
+
+ return bd, xmlroot
+
+ def get_key(self, key1, key2=False):
+ """
+ Get key from the BladedData CDATA section and format to int32 or
+ float32 numpy arrays if possible withouth precision loss.
+
+ Parameters
+ ----------
+ key1 : str
+ DESCRIPTION.
+ key2 : str, optional
+ DESCRIPTION. The default is False.
+
+
+ Returns
+ -------
+ numpy.array
+ Values from key1/key2 formatted as a numpy array. Converted to
+ numpy.int32, numpy.float32 if possible withouth precision loss,
+ otherwise an object array is returned.
+
+ """
+
+ if key1 not in self.bd:
+ raise KeyError(f'{key1} not found in BLADED file.')
+
+ if key2 is not False:
+ if key2 not in self.bd[key1]:
+ raise KeyError(f'{key2} not found in MSTART {key1} of BLADED file.')
+ data = self.bd[key1][key2]
+ else:
+ data = self.bd[key1]
+
+ # in case we defined a mstart key
+ if isinstance(data, dict):
+ return data
+
+ # i ,j = len(data), len(data[0])
+
+ # FIXME: this is a very expensive way of converting it, but it might
+ # not matter at all since very little model data is actually considered
+ data_arr = np.array(data)
+ try:
+ data_arr = data_arr.astype(np.int32)
+ except ValueError:
+ try:
+ data_arr = data_arr.astype(np.float32)
+ except ValueError:
+ pass
+
+ return data_arr
+
+ # return np.array(data, dtype=np.float32)
+
+ # if isinstance(data[0], list) and len(data[0]) == 1:
+ # data = float(data[0])
+ # if isinstance(data[0], list) and len(data[0]) > 1:
+ # data_arr = np.array(data, dtype=np.float32)
+
+ def set_attr_and_check(self):
+ """Check that BGEOMMB, BMASSMB, BSTIFFMB has indeed the same node
+ repated twice every time.
+ """
+
+ # self.bd['BGEOMMB'].keys(), but only those relevant
+ keysg = ['RJ', 'DIST', 'REF_X', 'REF_Y', 'CHORD', 'TWIST', 'CE_X',
+ 'CE_Y', 'BTHICK', 'FOIL', 'MOVING']
+ nbe = self.get_key('BGEOMMB', 'NBE')[0,0]
+
+ # self.bd['BMASSMB'].keys(), but only those relevant
+ keysm = ['CM_X', 'CM_Y', 'MASS', 'SINER', 'RGRATIO', 'BETA_M']
+
+ # self.bd['BSTIFFMB'].keys(), but only those relevant
+ keyss = ['EIFLAP', 'EIEDGE', 'BETA_S', 'GJ', 'CS_X', 'CS_Y', 'GAFLAP',
+ 'GAEDGE']
+
+ mkeys = ['BGEOMMB', 'BMASSMB', 'BSTIFFMB']
+ for mkey, keys in zip(mkeys, [keysg, keysm, keyss]):
+ for key in keys:
+ res = self.get_key(mkey, key)
+ try:
+ assert np.allclose(res[0,0::2], res[0,1::2])
+ except TypeError:
+ # allclose doesn't make sense for text arrays
+ assert np.compare_chararrays(res[0,0::2], res[0,1::2],
+ '==', True).all()
+ assert res.shape[1]==nbe
+ if hasattr(self, key.lower()):
+ raise UserWarning(key, 'already exists')
+ setattr(self, key.lower(), res[0,0::2])
+
+ def print_xml_tree(self, fname):
+ """For discovery purposes: print full tree + values/text
+
+
+ Parameters
+ ----------
+ fname : TYPE
+ DESCRIPTION.
+
+ Returns
+ -------
+ None.
+
+ """
+
+ # def print_rec(root):
+ # # if hasattr(root, 'getparent'):
+ # # print(root.getparent().tag.title, end='.')
+ # # print_rec(root.getparent())
+ # for el in root.getchildren():
+ # print(print_rec(el))
+
+ # def print_root_tree(root):
+ # root.getparent()
+ # print()
+
+ # tree = etree.fromstring(xml_str)
+ # els = tree.xpath('/')
+ # for el in els:
+ # print(el)
+
+ # tree = etree.fromstring(xml_str)
+
+ # xmlroot = objectify.fromstring(xml_str)
+ # with open(fname+'.structure', 'w') as f:
+ # for line in xmlroot.descendantpaths():
+ # f.write(line + '\n')
+
+ # Recursive XML parsing python using ElementTree
+ # https://stackoverflow.com/q/28194703/3156685
+
+ roottree = self.xmlroot.getroottree()
+ def print_tree_recursive(root):
+ # print(roottree.getpath(root), end=' : ')
+ # print(root.tag.title())
+ f.write(roottree.getpath(root) + ' : ')
+ f.writelines(root.values())
+ if root.text is not None:
+ f.write(' : ' + root.text)
+ f.write('\n')
+ for elem in root.getchildren():
+ print_tree_recursive(elem)
+ with open(fname+'.structure.value', 'w') as f:
+ for el in self.xmlroot.getchildren():
+ if el.tag.title()!='Bladeddata':
+ print_tree_recursive(el)
+
+ # def get_frequencies(self):
+ # """
+ # """
+
+ # blades = self.xmlroot.NewGUI.Turbine.Blades.FlexibilityModel
+
+ # blades.Settings.ModesWithDampingDefined
+ # blades.PartModes # with lots of BladeModeContainer
+ # blades.WholeBladeModes
+ # blades.WholeBladeModes.WholeBladeMode
+
+ # blades.WholeBladeModes.WholeBladeMode.Name
+ # blades.WholeBladeModes.WholeBladeMode.Frequency
+ # blades.WholeBladeModes.WholeBladeMode.Damping
+ # blades.WholeBladeModes.WholeBladeMode.Components
+ # len(blades.WholeBladeModes.WholeBladeMode.Components.getchildren())
+ # # 60 elements
+
+ # for wholeblademode in blades.WholeBladeModes.iterchildren():
+ # print(wholeblademode.Name)
+ # print(wholeblademode.Frequency, wholeblademode.Damping)
+
+ # tower = self.xmlroot.NewGUI.Turbine.SupportStructure.FlexibilityModel
+ # for towermode in tower.Modes.getchildren():
+ # print(towermode.Description)
+ # towermode.Frequency
+ # towermode.Damping
diff --git a/wetb/bladed/template.htc b/wetb/bladed/template.htc
new file mode 100644
index 0000000000000000000000000000000000000000..3d58aab3345e23bef7a58a94ed8da02e95d0574e
--- /dev/null
+++ b/wetb/bladed/template.htc
@@ -0,0 +1,324 @@
+;
+begin simulation;
+ time_stop 1;
+ solvertype 1 ; (newmark)
+ on_no_convergence continue ;
+ convergence_limits 1E3 1.0 1E-7 ;
+ logfile case0.log ;
+ begin newmark;
+ deltat 0.02;
+ end newmark;
+end simulation;
+;----------------------------------------------------------------------------------------------------------------------
+begin new_htc_structure;
+;----------------------------------------------------------------------------------------------------------------------
+ begin main_body;
+ name tower ;
+ type timoschenko ;
+ nbodies 1 ;
+ node_distribution c2_def ;
+ damping_posdef 0.0 0.0 0.0 0.0E-00 0.0E-00 0.0E-00;
+ begin timoschenko_input;
+ filename template.st;
+ set 1 1 ;
+ end timoschenko_input;
+ begin c2_def;
+ nsec 2;
+ sec 1 0 0 0.00 0 ;
+ sec 2 0 0 -100.00 0 ;
+ end c2_def ;
+ end main_body;
+;
+ begin main_body;
+ name towertop ;
+ type timoschenko ;
+ nbodies 1 ;
+ node_distribution c2_def ;
+ damping_posdef 0.0 0.0 0.0 0.0E-00 0.0E-00 0.0E-00;
+ begin timoschenko_input;
+ filename template.st;
+ set 1 1 ;
+ end timoschenko_input;
+ begin c2_def;
+ nsec 2;
+ sec 1 0.0 0.0 0.0 0.0 ;
+ sec 2 0.0 0.0 -1.00 0.0 ;
+ end c2_def ;
+ end main_body;
+;
+ begin main_body;
+ name shaft ;
+ type timoschenko ;
+ nbodies 1 ;
+ node_distribution c2_def ;
+ damping_posdef 0.0 0.0 0.0 0.0E-00 0.0E-00 0.0E-00;
+ begin timoschenko_input;
+ filename template.st;
+ set 1 1 ;
+ end timoschenko_input;
+ begin c2_def;
+ nsec 2;
+ sec 1 0.0 0.0 0.0 0.0 ;
+ sec 2 0.0 0.0 1.0 0.0 ;
+ end c2_def ;
+ end main_body;
+;
+ begin main_body;
+ name hub1 ;
+ type timoschenko ;
+ nbodies 1 ;
+ node_distribution c2_def ;
+ damping_posdef 0.0 0.0 0.0 0.0E-00 0.0E-00 0.0E-00;
+ begin timoschenko_input;
+ filename template.st;
+ set 1 1 ;
+ end timoschenko_input;
+ begin c2_def;
+ nsec 2;
+ sec 1 0.0 0.0 0.0 0.0 ;
+ sec 2 0.0 0.0 1.0 0.0 ;
+ end c2_def ;
+ end main_body;
+;
+ begin main_body;
+ name hub2 ;
+ copy_main_body hub1;
+ end main_body;
+;
+ begin main_body;
+ name hub3 ;
+ copy_main_body hub1 ;
+ end main_body;
+;
+ begin main_body;
+ name blade1 ;
+ type timoschenko ;
+ nbodies 1 ;
+ node_distribution c2_def;
+ damping_posdef 0.0 0.0 0.0 0.0E-00 0.0E-00 0.0E-00;
+ begin timoschenko_input ;
+ filename template.st;
+ set 1 1 ;
+ end timoschenko_input;
+ begin c2_def;
+ nsec 2;
+ sec 1 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 ;
+ sec 2 0.00000E+00 0.00000E+00 8.00000E+01 0.00000E+00 ;
+ end c2_def ;
+ end main_body;
+;
+ begin main_body;
+ name blade2 ;
+ copy_main_body blade1;
+ end main_body;
+;
+ begin main_body;
+ name blade3 ;
+ copy_main_body blade1 ;
+ end main_body;
+;----------------------------------------------------------------------------------------------------------------------
+;
+ begin orientation;
+ begin base;
+ mbdy tower;
+ inipos 0.0 0.0 0.0 ;
+ body_eulerang 0.0 0.0 0.0;
+ end base;
+ begin relative;
+ mbdy1 tower last;
+ mbdy2 towertop 1;
+ mbdy2_eulerang 0.0 0.0 0.0;
+ end relative;
+ begin relative;
+ mbdy1 towertop last;
+ mbdy2 shaft 1;
+ mbdy2_eulerang 90.0 0.0 0.0;
+ mbdy2_eulerang 0.0 0.0 0.0; [tilt] 0.0 0.0
+ mbdy2_eulerang 0.0 0.0 0.0; 0.0 0.0 [azi]
+ mbdy2_ini_rotvec_d1 0.0 0.0 -1.0 0.0 ; [init_wr];
+ end relative;
+ begin relative;
+ mbdy1 shaft last;
+ mbdy2 hub1 1;
+ mbdy2_eulerang -90.0 0.0 0.0;
+ mbdy2_eulerang 0.0 180.0 0.0;
+ mbdy2_eulerang 2.5 0.0 0.0; [cone] 0.0 0.0
+ end relative;
+ begin relative;
+ mbdy1 shaft last;
+ mbdy2 hub2 1;
+ mbdy2_eulerang -90.0 0.0 0.0;
+ mbdy2_eulerang 0.0 60.0 0.0;
+ mbdy2_eulerang 2.5 0.0 0.0;
+ end relative;
+ begin relative;
+ mbdy1 shaft last;
+ mbdy2 hub3 1;
+ mbdy2_eulerang -90.0 0.0 0.0;
+ mbdy2_eulerang 0.0 -60.0 0.0;
+ mbdy2_eulerang 2.5 0.0 0.0;
+ end relative;
+ begin relative;
+ mbdy1 hub1 last;
+ mbdy2 blade1 1;
+ mbdy2_eulerang 0.0 0.0 0;
+ end relative;
+ begin relative;
+ mbdy1 hub2 last;
+ mbdy2 blade2 1;
+ mbdy2_eulerang 0.0 0.0 0.0;
+ end relative;
+ begin relative;
+ mbdy1 hub3 last;
+ mbdy2 blade3 1;
+ mbdy2_eulerang 0.0 0.0 0.0;
+ end relative;
+ end orientation;
+;----------------------------------------------------------------------------------------------------------------------
+ begin constraint;
+ begin fix0;
+ body tower;
+ end fix0;
+ begin fix1;
+ mbdy1 tower last ;
+ mbdy2 towertop 1;
+ end fix1;
+ begin bearing1;
+ name shaft_rot;
+ mbdy1 towertop last;
+ mbdy2 shaft 1;
+ bearing_vector 2 0.0 0.0 -1.0; (0=global.1=mbdy1.2=mbdy2)
+ end bearing1;
+; begin bearing3;
+; name shaft_rot;
+; mbdy1 towertop last;
+; mbdy2 shaft 1;
+; bearing_vector 2 0.0 0.0 -1.0; (0=global.1=mbdy1.2=mbdy2)
+; omegas 0.0 ;
+; end bearing3;
+ begin fix1;
+ mbdy1 shaft last ;
+ mbdy2 hub1 1;
+ end fix1;
+ begin fix1;
+ mbdy1 shaft last ;
+ mbdy2 hub2 1;
+ end fix1;
+ begin fix1;
+ mbdy1 shaft last ;
+ mbdy2 hub3 1;
+ end fix1;
+ begin bearing2;
+ name pitch1;
+ mbdy1 hub1 last;
+ mbdy2 blade1 1;
+ bearing_vector 2 0.0 0.0 -1.0;
+ end bearing2;
+ begin bearing2;
+ name pitch2;
+ mbdy1 hub2 last;
+ mbdy2 blade2 1;
+ bearing_vector 2 0.0 0.0 -1.0;
+ end bearing2;
+ begin bearing2;
+ name pitch3;
+ mbdy1 hub3 last;
+ mbdy2 blade3 1;
+ bearing_vector 2 0.0 0.0 -1.0;
+ end bearing2;
+ end constraint;
+end new_htc_structure;
+;----------------------------------------------------------------------------------------------------------------------
+begin wind ;
+ density 1.225 ;
+ wsp 1;
+ tint 0.01;
+ horizontal_input 1 ;
+ windfield_rotations 0.0 0.0 0.0 ; yaw, tilt, rotation
+ center_pos0 0.0 0.0 -119 ; hub heigth
+ shear_format 3 0.2;
+ turb_format 0 ; 0=none, 1=mann,2=flex
+ tower_shadow_method 3 ; 0=none, 1=potential flow, 2=jet
+; scale_time_start 100;
+; wind_ramp_factor 0.0 50 0.0 1.0 ;
+;
+ begin tower_shadow_potential_2;
+ tower_mbdy_link tower;
+ nsec 2;
+ radius 0.0 4.15 ;
+ radius 115.63 2.75 ;
+ end tower_shadow_potential_2;
+end wind;
+;
+begin aerodrag ;
+ begin aerodrag_element ;
+ mbdy_name tower;
+ aerodrag_sections uniform 10 ;
+ nsec 2 ;
+ sec 0.0 0.6 1.0 ; tower bottom
+ sec 100.0 0.6 1.0 ; tower top
+ end aerodrag_element;
+;
+ begin aerodrag_element ; Nacelle drag side
+ mbdy_name shaft;
+ aerodrag_sections uniform 2 ;
+ nsec 2 ;
+ sec 0.0 0.8 10.0 ;
+ sec 1.00 0.8 10.0 ;
+ end aerodrag_element;
+end aerodrag;
+;
+begin aero ;
+ nblades 3;
+ hub_vec shaft -3 ;
+ link 1 mbdy_c2_def blade1;
+ link 2 mbdy_c2_def blade2;
+ link 3 mbdy_c2_def blade3;
+ ae_filename template.ae ;
+ pc_filename template.pc ;
+ induction_method 0 ; 0=none, 1=normal
+ aerocalc_method 1 ; 0=ingen aerodynamic, 1=med aerodynamic
+ aerosections 2;
+ ae_sets 1 1 1;
+ tiploss_method 1 ; 0=none, 1=prandtl
+ dynstall_method 2 ; 0=none, 1=stig øye method,2=mhh method
+end aero ;
+;----------------------------------------------------------------------------------------------------------------------
+begin output;
+ filename case0;
+ time 0 100;
+ data_format gtsdf;
+ buffer 9999;
+;
+ general time;
+ constraint bearing1 shaft_rot 2; angle and angle velocity
+ constraint bearing2 pitch1 5; angle and angle velocity
+ constraint bearing2 pitch2 5; angle and angle velocity
+ constraint bearing2 pitch3 5; angle and angle velocity
+ aero omega ;
+ aero torque;
+ aero power;
+ aero thrust;
+ wind free_wind 1 0.0 0.0 -100;
+ ; Moments:
+ mbdy momentvec tower 1 1 tower # tower base ;
+ mbdy momentvec tower 1 2 tower # tower yaw bearing ;
+ mbdy momentvec shaft 1 1 shaft # main bearing ;
+ mbdy momentvec blade1 1 1 blade1 # blade 1 root ;
+ mbdy momentvec blade2 1 1 blade2 # blade 2 root ;
+ mbdy momentvec blade3 1 1 blade3 # blade 3 root ;
+ mbdy momentvec blade1 1 2 local # blade 1 50% local e coo ;
+ mbdy momentvec blade2 1 2 local # blade 2 50% local e coo ;
+ mbdy momentvec blade3 1 2 local # blade 3 50% local e coo ;
+ ; Displacements and accellerations
+ mbdy state pos tower 1 1.0 global only 1 # Tower top FA displ;
+ mbdy state pos tower 1 1.0 global only 2 # Tower top SS displ;
+ mbdy state acc tower 1 1.0 global only 1 # Tower top FA acc;
+ mbdy state acc tower 1 1.0 global only 2 # Tower top SS acc;
+ mbdy state pos blade1 1 1.0 global # gl blade 1 tip pos ;
+ mbdy state pos blade2 1 1.0 global # gl blade 2 tip pos ;
+ mbdy state pos blade3 1 1.0 global # gl blade 3 tip pos ;
+ mbdy state pos blade1 1 1.0 blade1 # blade 1 tip pos ;
+end output;
+;
+exit;
diff --git a/wetb/bladed/template.st b/wetb/bladed/template.st
new file mode 100644
index 0000000000000000000000000000000000000000..a4105bad9aaeec43c845ca8cb144e622d1fe0842
--- /dev/null
+++ b/wetb/bladed/template.st
@@ -0,0 +1,10 @@
+1 number of sets, Nset
+-----------------
+#1 Nset number 1
+================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================
+ r [0] m [1] x_cg [2] y_cg [3] ri_x [4] ri_y [5] x_sh [6] y_sh [7] E [8] G [9] I_x [10] I_y [11] K [12] k_x [13] k_y [14] A [15] pitch [16] x_e [17] y_e [18]
+================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================
+$1 2 subset number 1
+ 0.000000000000000e+00 0.500000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.500000000000000e+00 0.500000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+14 1.000000000000000e+13 1.000000000000000e+00 1.000000000000000e+00 1.000000000000000e+00 1.000000000000000e+00 1.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
+ 1.000000000000000e+00 0.500000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.500000000000000e+00 0.500000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+14 1.000000000000000e+13 1.000000000000000e+00 1.000000000000000e+00 1.000000000000000e+00 1.000000000000000e+00 1.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
+