Simulations.py

        load case. Starting from Mx and My moments, the other cross-sectional
        forces are identified.

        Parameters
        ----------

        sig : list, time-series signal

        ch_list : list, list of channels for enevelope computation

        int_env : boolean, default=False
            If the logic parameter is True, the function will interpolate the
            envelope on a given number of points

        Nx : int, default=300
            Number of points for the envelope interpolation

        Returns
        -------

        envelope : dictionary,
            The dictionary has entries refered to the channels selected.
            Inside the dictonary under each entry there is a matrix with 6
            columns, each for the sectional forces and moments

        """

        envelope= {}
        for ch in ch_list:
            chi0 = self.res.ch_dict[ch[0]]['chi']
            chi1 = self.res.ch_dict[ch[1]]['chi']
            s0 = np.array(sig[:, chi0]).reshape(-1, 1)
            s1 = np.array(sig[:, chi1]).reshape(-1, 1)
            # Compute a Convex Hull, the vertices number varies according to
            # the shape of the poligon
            cloud =  np.append(s0, s1, axis=1)
            hull = scipy.spatial.ConvexHull(cloud)
            closed_contour = np.append(cloud[hull.vertices,:],
                                       cloud[hull.vertices[0],:].reshape(1,2),
                                       axis=0)

            # Interpolate envelope for a given number of points
            if int_env:
                _,_,_,closed_contour_int = int_envelope(closed_contour[:,0],
                                                        closed_contour[:,1],Nx)


            # Based on Mx and My envelope, the other cross-sectional moments
            # and forces components are identified and appended to the initial
            # envelope
            for ich in range(2, len(ch)):
                chix = self.res.ch_dict[ch[ich]]['chi']
                s0 = np.array(sig[hull.vertices, chix]).reshape(-1, 1)
                s1 = np.array(sig[hull.vertices[0], chix]).reshape(-1, 1)
                s0 = np.append(s0, s1, axis=0)
                closed_contour = np.append(closed_contour, s0, axis=1)
                if int_env:
                    _,_,_,extra_sensor = int_envelope(closed_contour[:,0],
                                                       closed_contour[:,ich],Nx)
                    es = np.atleast_2d(np.array(extra_sensor[:,1])).T
                    closed_contour_int = np.append(closed_contour_int,es,axis=1)

            if int_env:
                envelope[ch[0]] = closed_contour_int
            else:
                envelope[ch[0]] = closed_contour

        return envelope

    def int_envelope(ch1,ch2,Nx):
        # Function to interpolate envelopes and output arrays of same length

        # Number of points is defined by Nx + 1, where the + 1 is needed to
        # close the curve

        upper = []
        lower = []

        indmax = np.argmax(ch1)
        indmin = np.argmin(ch1)
        if indmax > indmin:
            lower = np.array([ch1[indmin:indmax+1],ch2[indmin:indmax+1]]).T
            upper = np.concatenate((np.array([ch1[indmax:],ch2[indmax:]]).T,\
                            np.array([ch1[:indmin+1],ch2[:indmin+1]]).T),axis=0)
        else:
            upper = np.array([ch1[indmax:indmin+1,:],ch2[indmax:indmin+1,:]]).T
            lower = np.concatenate((np.array([ch1[indmin:],ch2[indmin:]]).T,\
                                np.array([ch1[:indmax+1],ch2[:indmax+1]]).T),axis=0)


        int_1 = np.linspace(min(min(upper[:,0]),min(lower[:,0])),\
                            max(max(upper[:,0]),max(upper[:,0])),Nx/2+1)
        upper = np.flipud(upper)
        int_2_up = np.interp(int_1,np.array(upper[:,0]),np.array(upper[:,1]))
        int_2_low = np.interp(int_1,np.array(lower[:,0]),np.array(lower[:,1]))

        int_env = np.concatenate((np.array([int_1[:-1],int_2_up[:-1]]).T,\
                                np.array([int_1[::-1],int_2_low[::-1]]).T),axis=0)

        return int_env

    def envelope(self, silent=False, ch_list=[], append=''):
        """
        Calculate envelopes and save them in a table.

        Parameters
        ----------


        Returns
        -------


        """
        # get some basic parameters required to calculate statistics
        try:
            case = list(self.cases.keys())[0]
        except IndexError:
            print('no cases to select so no statistics, aborting ...')
            return None

        post_dir = self.cases[case]['[post_dir]']
        sim_id = self.cases[case]['[sim_id]']

        if not silent:
            nrcases = len(self.cases)
            print('='*79)
            print('statistics for %s, nr cases: %i' % (sim_id, nrcases))

        fname = os.path.join(post_dir, sim_id + '_envelope' + append + '.h5')
        h5f = tbl.openFile(fname, mode="w", title=str(sim_id),
                           filters=tbl.Filters(complevel=9))

        # Create a new group under "/" (root)
        for ii, (cname, case) in enumerate(self.cases.items()):

            groupname = str(cname[:-4])
            groupname = groupname.replace('-', '_')
            h5f.createGroup("/", groupname)
            ctab = getattr(h5f.root, groupname)

            if not silent:
                pc = '%6.2f' % (float(ii)*100.0/float(nrcases))
                pc += ' %'
                print('envelope progress: %4i/%i %s' % (ii, nrcases, pc))

            self.load_result_file(case)

            envelope = self.compute_envelope(self.sig, ch_list)

            for ch_id in ch_list:
                h5f.createTable(ctab, str(ch_id[0].replace('-', '_')),
                                EnvelopeClass.section,
                                title=str(ch_id[0].replace('-', '_')))
                csv_table = getattr(ctab, str(ch_id[0].replace('-', '_')))
                tablerow = csv_table.row
                for row in envelope[ch_id[0]]:
                    tablerow['Mx'] = float(row[0])
                    tablerow['My'] = float(row[1])
                    if len(row)>2:
                        tablerow['Mz'] = float(row[2])
                        if len(row)>3:
                            tablerow['Fx'] = float(row[3])
                            tablerow['Fy'] = float(row[4])
                            tablerow['Fz'] = float(row[5])
                        else:
                            tablerow['Fx'] = 0.0
                            tablerow['Fy'] = 0.0
                            tablerow['Fz'] = 0.0
                    else:
                        tablerow['Mz'] = 0.0
                        tablerow['Fx'] = 0.0
                        tablerow['Fy'] = 0.0
                        tablerow['Fz'] = 0.0
                    tablerow.append()
                csv_table.flush()
        h5f.close()


class EnvelopeClass(object):
    """
    Class with the definition of the table for the envelope results
    """
    class section(tbl.IsDescription):

        Mx = tbl.Float32Col()
        My = tbl.Float32Col()
        Mz = tbl.Float32Col()
        Fx = tbl.Float32Col()
        Fy = tbl.Float32Col()
        Fz = tbl.Float32Col()


# TODO: implement this
class Results(object):
    """
    Move all Hawc2io to here? NO: this should be the wrapper, to interface
    the htc_dict with the io functions

    There should be a bare metal module/class for those who only want basic
    python support for HAWC2 result files and/or launching simulations.

    How to properly design this module? Change each class into a module? Or
    leave like this?
    """

    # OK, for now use this to do operations on HAWC2 results files

    def __init___(self):
        """
        """
        pass

    def m_equiv(self, st_arr, load, pos):
        r"""Centrifugal corrected equivalent moment

        Convert beam loading into a single equivalent bending moment. Note that
        this is dependent on the location in the cross section. Due to the
        way we measure the strain on the blade and how we did the calibration
        of those sensors.

        .. math::

            \epsilon = \frac{M_{x_{equiv}}y}{EI_{xx}} = \frac{M_x y}{EI_{xx}}
            + \frac{M_y x}{EI_{yy}} + \frac{F_z}{EA}

            M_{x_{equiv}} = M_x + \frac{I_{xx}}{I_{yy}} M_y \frac{x}{y}
            + \frac{I_{xx}}{Ay} F_z

        Parameters
        ----------

        st_arr : np.ndarray(19)
            Only one line of the st_arr is allowed and it should correspond
            to the correct radial position of the strain gauge.

        load : list(6)
            list containing the load time series of following components
            .. math:: load = F_x, F_y, F_z, M_x, M_y, M_z
            and where each component is an ndarray(m)

        pos : np.ndarray(2)
            x,y position wrt neutral axis in the cross section for which the
            equivalent load should be calculated

        Returns
        -------

        m_eq : ndarray(m)
            Equivalent load, see main title

        """

        F_z = load[2]
        M_x = load[3]
        M_y = load[4]

        x, y = pos[0], pos[1]

        A = st_arr[ModelData.st_headers.A]
        I_xx = st_arr[ModelData.st_headers.Ixx]
        I_yy = st_arr[ModelData.st_headers.Iyy]

        M_x_equiv = M_x + ( (I_xx/I_yy)*M_y*(x/y) ) + ( F_z*I_xx/(A*y) )
        # or ignore edgewise moment
        #M_x_equiv = M_x + ( F_z*I_xx/(A*y) )

        return M_x_equiv


class MannTurb64(prepost.PBSScript):
    """
    alfaeps, L, gamma, seed, nr_u, nr_v, nr_w, du, dv, dw high_freq_comp
    mann_turb_x64.exe fname 1.0 29.4 3.0 1209 256 32 32 2.0 5 5 true
    """

    def __init__(self, silent=False):
        super(MannTurb64, self).__init__()
        self.exe = 'time wine mann_turb_x64.exe'
        # PBS configuration
        self.umask = '0003'
        self.walltime = '00:59:59'
        self.queue = 'workq'
        self.lnodes = '1'
        self.ppn = '1'
        self.silent = silent
        self.pbs_in_dir = 'pbs_in_turb/'

    def gen_pbs(self, cases):

        case0 = cases[list(cases.keys())[0]]
        # make sure the path's end with a trailing separator, why??
        self.pbsworkdir = os.path.join(case0['[run_dir]'], '')
        if not self.silent:
            print('\nStart creating PBS files for turbulence with Mann64...')
        for cname, case in cases.items():

            # only relevant for cases with turbulence
            if '[tu_model]' in case and int(case['[tu_model]']) == 0:
                continue
            if '[Turb base name]' not in case:
                continue

            base_name = case['[Turb base name]']
            # pbs_in/out dir can contain subdirs, only take the inner directory
            out_base = misc.path_split_dirs(case['[pbs_out_dir]'])[0]
            turb = case['[turb_dir]']

            self.path_pbs_e = os.path.join(out_base, turb, base_name + '.err')
            self.path_pbs_o = os.path.join(out_base, turb, base_name + '.out')
            self.path_pbs_i = os.path.join(self.pbs_in_dir, base_name + '.p')

            if case['[turb_db_dir]'] is not None:
                self.prelude = 'cd %s' % case['[turb_db_dir]']
            else:
                self.prelude = 'cd %s' % case['[turb_dir]']

            # alfaeps, L, gamma, seed, nr_u, nr_v, nr_w, du, dv, dw high_freq_comp
            rpl = (float(case['[MannAlfaEpsilon]']),
                   float(case['[MannL]']),
                   float(case['[MannGamma]']),
                   int(case['[tu_seed]']),
                   int(case['[turb_nr_u]']),
                   int(case['[turb_nr_v]']),
                   int(case['[turb_nr_w]']),
                   float(case['[turb_dx]']),
                   float(case['[turb_dy]']),
                   float(case['[turb_dz]']),
                   int(case['[high_freq_comp]']))
            params = '%1.6f %1.6f %1.6f %i %i %i %i %1.4f %1.4f %1.4f %i' % rpl
            self.execution = '%s %s %s' % (self.exe, base_name, params)
            self.create(check_dirs=True)


def eigenbody(cases, debug=False):
    """
    Read HAWC2 body eigenalysis result file
    =======================================

    This is basically a cases convience wrapper around Hawc2io.ReadEigenBody

    Parameters
    ----------

    cases : dict{ case : dict{tag : value} }
        Dictionary where each case is a key and its value a dictionary
        holding all the tags/value pairs as used for that case.

    Returns
    -------

    cases : dict{ case : dict{tag : value} }
        Dictionary where each case is a key and its value a dictionary
        holding all the tags/value pairs as used for that case. For each
        case, it is updated with the results, results2 of the eigenvalue
        analysis performed for each body using the following respective
        tags: [eigen_body_results] and [eigen_body_results2].

    """

    #Body data for body number : 3 with the name :nacelle
    #Results:         fd [Hz]       fn [Hz]       log.decr [%]
    #Mode nr:  1:   1.45388E-21    1.74896E-03    6.28319E+02

    for case in cases:
        # tags for the current case
        tags = cases[case]
        file_path = os.path.join(tags['[run_dir]'], tags['[eigenfreq_dir]'])
        # FIXME: do not assuem anything about the file name here, should be
        # fully defined in the tags/dataframe
        file_name = tags['[case_id]'] + '_body_eigen'
        # and load the eigenfrequency body results
        results, results2 = windIO.ReadEigenBody(file_path, file_name,
                                                  nrmodes=10)
        # add them to the htc_dict
        cases[case]['[eigen_body_results]'] = results
        cases[case]['[eigen_body_results2]'] = results2

    return cases

def eigenstructure(cases, debug=False):
    """
    Read HAWC2 structure eigenalysis result file
    ============================================

    This is basically a cases convience wrapper around
    windIO.ReadEigenStructure

    Parameters
    ----------

    cases : dict{ case : dict{tag : value} }
        Dictionary where each case is a key and its value a dictionary
        holding all the tags/value pairs as used for that case.

    Returns
    -------

    cases : dict{ case : dict{tag : value} }
        Dictionary where each case is a key and its value a dictionary
        holding all the tags/value pairs as used for that case. For each
        case, it is updated with the modes_arr of the eigenvalue
        analysis performed for the structure.
        The modes array (ndarray(3,n)) holds fd, fn and damping.
    """

    for case in cases:
        # tags for the current case
        tags = cases[case]
        file_path = os.path.join(tags['[run_dir]'], tags['[eigenfreq_dir]'])
        # FIXME: do not assuem anything about the file name here, should be
        # fully defined in the tags/dataframe
        file_name = tags['[case_id]'] + '_strc_eigen'
        # and load the eigenfrequency structure results
        modes = windIO.ReadEigenStructure(file_path, file_name, max_modes=500)
        # add them to the htc_dict
        cases[case]['[eigen_structure]'] = modes

    return cases

if __name__ == '__main__':
    pass