#!/usr/bin/env python
# encoding: utf-8
"""
towerstruc.py

Created by Andrew Ning on 2012-01-20.
Copyright (c) NREL. All rights reserved.

HISTORY:  2012 created
          -7/2014:  R.D. Bugs found in the call to shellBucklingEurocode from towerwithFrame3DD. Fixed.
                    Also set_as_top added.
          -10/2014: R.D. Merged back with some changes Andrew did on his end.
          -12/2014: A.N. fixed some errors from the merge (redundant drag calc).  pep8 compliance.  removed several unneccesary variables and imports (including set_as_top)
          - 6/2015: A.N. major rewrite.  removed pBEAM.  can add spring stiffness anywhere.  can add mass anywhere.
            can use different material props throughout.
          - 7/2015 : R.D. modified to use commonse modules.
          - 1/2018 : G.B. modified for easier use with other modules, reducing user input burden, and shifting more to commonse
 """

from __future__ import print_function

import numpy as np
from openmdao.api import ExplicitComponent, Group, Problem, IndepVarComp

from wisdem.commonse.WindWaveDrag import AeroHydroLoads, CylinderWindDrag, CylinderWaveDrag

from wisdem.commonse.environment import WindBase, WaveBase, LinearWaves, TowerSoil, PowerWind, LogWind
from wisdem.commonse.tube import CylindricalShellProperties
from wisdem.commonse.utilities import assembleI, unassembleI, nodal2sectional, interp_with_deriv
from wisdem.commonse import gravity, eps, NFREQ

from wisdem.commonse.vertical_cylinder import CylinderDiscretization, CylinderMass, CylinderFrame3DD

import wisdem.commonse.UtilizationSupplement as Util


def find_nearest(array,value):
    return (np.abs(array-value)).argmin() 


# -----------------
#  Components
# -----------------

class MonopileFoundation(ExplicitComponent):
    def initialize(self):
        self.options.declare('nPoints')
        self.options.declare('monopile')
    
    def setup(self):
        nPoints = self.options['nPoints']

        self.add_input('tower_section_height', np.zeros(nPoints-1), units='m', desc='parameterized section heights along cylinder')
        self.add_input('tower_outer_diameter', np.zeros(nPoints), units='m', desc='cylinder diameter at corresponding locations')
        self.add_input('tower_wall_thickness', np.zeros(nPoints-1), units='m', desc='shell thickness at corresponding locations')
        self.add_input('suctionpile_depth', 0.0, units='m', desc='depth of foundation in the soil')
        self.add_input('suctionpile_depth_diam_ratio', 0.0, desc='ratio of sunction pile depth to mudline monopile diameter')
        self.add_input('foundation_height', 0.0, units='m', desc='height of foundation (0.0 for land, -water_depth for fixed bottom)')

        nadd = 1 if self.options['monopile'] else 0
        self.add_output('section_height_out', np.zeros(nPoints-1+nadd), units='m', desc='parameterized section heights along cylinder')
        self.add_output('outer_diameter_out', np.zeros(nPoints+nadd), units='m', desc='cylinder diameter at corresponding locations')
        self.add_output('wall_thickness_out', np.zeros(nPoints-1+nadd), units='m', desc='shell thickness at corresponding locations')
        self.add_output('foundation_height_out', 0.0, units='m', desc='height of suction pile bottom (0.0 for land, -water_depth-pile for fixed bottom)')

    def compute(self, inputs, outputs):
        if self.options['monopile']:
            # Determine suction pile depth
            pile = inputs['suctionpile_depth']
            if pile == 0.0:
                pile = inputs['suctionpile_depth_diam_ratio'] * inputs['tower_outer_diameter'][0]

            # Gravity foundations will have 0 suction depth, but we will still add an extra point for consistency
            pile = np.maximum(0.1, pile)
            
            outputs['section_height_out'] = np.r_[pile, inputs['tower_section_height']]
            outputs['outer_diameter_out'] = np.r_[inputs['tower_outer_diameter'][0], inputs['tower_outer_diameter']]
            outputs['wall_thickness_out'] = np.r_[inputs['tower_wall_thickness'][0], inputs['tower_wall_thickness']]
            outputs['foundation_height_out'] = inputs['foundation_height'] - pile
        else:
            outputs['section_height_out'] = inputs['tower_section_height']
            outputs['outer_diameter_out'] = inputs['tower_outer_diameter']
            outputs['wall_thickness_out'] = inputs['tower_wall_thickness']
            outputs['foundation_height_out'] = inputs['foundation_height']
            

class TowerDiscretization(ExplicitComponent):
    def initialize(self):
        self.options.declare('nPoints')
    
    def setup(self):
        nPoints = self.options['nPoints']

        self.add_input('hub_height', val=0.0, units='m', desc='diameter at tower base')
        self.add_input('z_param', np.zeros(nPoints), units='m', desc='parameterized locations along tower, linear lofting between')

        self.add_output('height_constraint', val=0.0, units='m', desc='mismatch between tower height and desired hub_height')

        self.declare_partials('height_constraint', ['hub_height','z_param'])
        
    def compute(self, inputs, outputs):
        outputs['height_constraint'] = inputs['hub_height'] - inputs['z_param'][-1]

    def compute_partials(self, inputs, J):
        nPoints = self.options['nPoints']
        
        J['height_constraint','hub_height'] = 1.
        J['height_constraint','z_param'] = np.zeros(nPoints)
        J['height_constraint','z_param'][-1] = -1.
        
        
        
class TowerMass(ExplicitComponent):

    def initialize(self):
        self.options.declare('nFull')
    
    def setup(self):
        nFull = self.options['nFull']
        
        self.add_input('cylinder_mass', val=np.zeros(nFull-1), units='kg', desc='Total cylinder mass')
        self.add_input('cylinder_cost', val=0.0, units='USD', desc='Total cylinder cost')
        self.add_input('cylinder_center_of_mass', val=0.0, units='m', desc='z-position of center of mass of cylinder')
        self.add_input('cylinder_section_center_of_mass', val=np.zeros(nFull-1), units='m', desc='z position of center of mass of each can in the cylinder')
        self.add_input('cylinder_I_base', np.zeros((6,)), units='kg*m**2', desc='mass moment of inertia of cylinder about base [xx yy zz xy xz yz]')

        self.add_input('transition_piece_height', 0.0, units='m', desc='point mass height of transition piece above water line')
        self.add_input('transition_piece_mass', 0.0, units='kg', desc='point mass of transition piece')
        self.add_input('gravity_foundation_mass', 0.0, units='kg', desc='extra mass of gravity foundation')
        self.add_input('foundation_height', 0.0, units='m', desc='height of foundation (0.0 for land, -water_depth for fixed bottom)')
        self.add_input('z_full', np.zeros(nFull), units='m', desc='parameterized locations along tower, linear lofting between')
        
        self.add_output('tower_raw_cost', val=0.0, units='USD', desc='Total tower cost')
        self.add_output('tower_mass', val=0.0, units='kg', desc='Total tower mass')
        self.add_output('tower_center_of_mass', val=0.0, units='m', desc='z-position of center of mass of tower')
        self.add_output('tower_section_center_of_mass', val=np.zeros(nFull-1), units='m', desc='z position of center of mass of each can in the tower')
        self.add_output('tower_I_base', np.zeros((6,)), units='kg*m**2', desc='mass moment of inertia of tower about base [xx yy zz xy xz yz]')

        self.add_output('monopile_mass', val=0.0, units='kg', desc='Mass of monopile from bottom of suction pile through transition piece')
        self.add_output('monopile_cost', val=0.0, units='USD', desc='Total monopile cost')
        self.add_output('monopile_length', val=0.0, units='m', desc='Length of monopile from bottom of suction pile through transition piece')
        
        self.declare_partials('tower_raw_cost', 'cylinder_cost')
        self.declare_partials('tower_mass', ['cylinder_mass','transition_piece_mass'])
        self.declare_partials('tower_center_of_mass', 'cylinder_center_of_mass')
        self.declare_partials('tower_section_center_of_mass', 'cylinder_section_center_of_mass')
        self.declare_partials('tower_I_base', 'cylinder_I_base')
        self.declare_partials('monopile_mass', ['cylinder_mass','z_full','transition_piece_height'])
        self.declare_partials('monopile_cost', ['cylinder_mass','z_full','transition_piece_height','cylinder_cost'])
        self.declare_partials('monopile_length', ['transition_piece_height','z_full'])

        self.J = {}
        
        
    def compute(self, inputs, outputs):
        outputs['tower_raw_cost']       = inputs['cylinder_cost']
        outputs['tower_mass']           = inputs['cylinder_mass'].sum()
        outputs['tower_center_of_mass'] = ( (inputs['cylinder_center_of_mass']*outputs['tower_mass'] +
                                             inputs['transition_piece_mass']*inputs['transition_piece_height'] +
                                             inputs['gravity_foundation_mass']*inputs['foundation_height']) /
                                            (outputs['tower_mass']+inputs['transition_piece_mass']+inputs['gravity_foundation_mass']) )
        outputs['tower_section_center_of_mass'] = inputs['cylinder_section_center_of_mass']
        outputs['tower_I_base']         = inputs['cylinder_I_base']

        outputs['monopile_mass'],dydx,dydxp,dydyp = interp_with_deriv(inputs['transition_piece_height'],
                                                                      inputs['z_full'],
                                                                      np.r_[0.0, np.cumsum(inputs['cylinder_mass'])])
        outputs['tower_mass']     -= outputs['monopile_mass']
        outputs['monopile_cost']   = inputs['cylinder_cost']*outputs['monopile_mass']/inputs['cylinder_mass'].sum()
        outputs['monopile_mass']  += inputs['transition_piece_mass'] + inputs['gravity_foundation_mass']
        outputs['monopile_length'] = inputs['transition_piece_height'] - inputs['z_full'][0]

        self.J = {}
        self.J['monopile_mass', 'z_full'] = dydxp[0,:]
        self.J['monopile_mass', 'cylinder_mass'] = dydyp[0,1:]
        self.J['monopile_mass', 'transition_piece_height'] = dydx[0,0]

        self.J['monopile_cost', 'z_full'] = inputs['cylinder_cost'] * self.J['monopile_mass', 'z_full'] / inputs['cylinder_mass'].sum()
        self.J['monopile_cost', 'cylinder_cost'] = outputs['monopile_mass']/ inputs['cylinder_mass']
        self.J['monopile_cost', 'cylinder_mass'] = inputs['cylinder_cost']*self.J['monopile_mass', 'cylinder_mass']/inputs['cylinder_mass'] - outputs['monopile_cost']/inputs['cylinder_mass']
        self.J['monopile_cost', 'transition_piece_height'] = inputs['cylinder_cost'] * self.J['monopile_mass', 'transition_piece_height'] / inputs['cylinder_mass']
        
        
    def compute_partials(self, inputs, J):
        J['tower_mass','cylinder_mass'] = np.ones(len(inputs['cylinder_mass'])) - self.J['monopile_mass', 'cylinder_mass']
        J['tower_mass','transition_piece_mass'] = 1.0
        J['tower_mass', 'z_full'] = -self.J['monopile_mass', 'z_full']
        J['tower_mass', 'transition_piece_height'] = -self.J['monopile_mass', 'transition_piece_height']

        J['tower_raw_cost','cylinder_cost'] = 1.0

        J['tower_center_of_mass','cylinder_center_of_mass'] = 1.0

        J['tower_section_center_of_mass','cylinder_section_center_of_mass'] = np.eye(len(inputs['cylinder_section_center_of_mass']))

        J['tower_I_base','cylinder_I_base'] = np.eye(len(inputs['cylinder_I_base']))

        J['monopile_mass', 'z_full'] = self.J['monopile_mass', 'z_full']
        J['monopile_mass', 'cylinder_mass'] = self.J['monopile_mass', 'cylinder_mass']
        J['monopile_mass', 'transition_piece_height'] = self.J['monopile_mass', 'transition_piece_height']
        J['monopile_mass', 'transition_piece_mass'] = 1.0
        J['monopile_mass', 'gravity_foundation_mass'] = 1.0

        J['monopile_cost', 'z_full'] = self.J['monopile_cost', 'z_full']
        J['monopile_cost', 'cylinder_cost'] = self.J['monopile_cost', 'cylinder_cost']
        J['monopile_cost', 'cylinder_mass'] = self.J['monopile_cost', 'cylinder_mass']
        J['monopile_cost', 'transition_piece_height'] = self.J['monopile_cost', 'transition_piece_height']

        J['monopile_length','transition_piece_height'] = 1.
        J['monopile_length','z_full'] = np.zeros(inputs['z_full'].size)
        J['monopile_length','z_full'][0] = -1.
        
        
        
class TurbineMass(ExplicitComponent):

    def setup(self):
        
        self.add_input('hub_height', val=0.0, units='m', desc='Hub-height')
        self.add_input('rna_mass', val=0.0, units='kg', desc='Total tower mass')
        self.add_input('rna_I', np.zeros((6,)), units='kg*m**2', desc='mass moment of inertia of rna about tower top [xx yy zz xy xz yz]')
        self.add_input('rna_cg', np.zeros((3,)), units='m', desc='xyz-location of rna cg relative to tower top')
        
        self.add_input('tower_mass', val=0.0, units='kg', desc='Total tower mass (not including monopile)')
        self.add_input('monopile_mass', val=0.0, units='kg', desc='Monopile mass')
        self.add_input('tower_center_of_mass', val=0.0, units='m', desc='z-position of center of mass of tower')
        self.add_input('tower_I_base', np.zeros((6,)), units='kg*m**2', desc='mass moment of inertia of tower about base [xx yy zz xy xz yz]')

        self.add_output('turbine_mass', val=0.0, units='kg', desc='Total mass of tower+rna')
        self.add_output('turbine_center_of_mass', val=np.zeros((3,)), units='m', desc='xyz-position of tower+rna center of mass')
        self.add_output('turbine_I_base', np.zeros((6,)), units='kg*m**2', desc='mass moment of inertia of tower about base [xx yy zz xy xz yz]')
       
        # Derivatives
        # self.declare_partials('*', '*', method='fd', form='central', step=1e-6)
        
    def compute(self, inputs, outputs):
        outputs['turbine_mass'] = inputs['rna_mass'] + inputs['tower_mass'] + inputs['monopile_mass']
        
        cg_rna   = inputs['rna_cg'] + np.array([0.0, 0.0, inputs['hub_height']])
        cg_tower = np.array([0.0, 0.0, inputs['tower_center_of_mass']])
        outputs['turbine_center_of_mass'] = (inputs['rna_mass']*cg_rna + inputs['tower_mass']*cg_tower) / outputs['turbine_mass']

        R       = cg_rna
        I_tower = assembleI(inputs['tower_I_base'])
        I_rna   = assembleI(inputs['rna_I']) + inputs['rna_mass']*(np.dot(R, R)*np.eye(3) - np.outer(R, R))
        outputs['turbine_I_base'] = unassembleI(I_tower + I_rna)

        

        

        
class TowerPreFrame(ExplicitComponent):
    def initialize(self):
        self.options.declare('nFull')
        self.options.declare('nPoints')
        self.options.declare('monopile', default=False)
    
    def setup(self):
        nPoints = self.options['nPoints']
        nFull   = self.options['nFull']
        nRefine = int( (nFull-1)/(nPoints-1) )
        
        self.add_input('z', np.zeros(nFull), units='m', desc='location along tower. start at bottom and go to top')
        self.add_input('d', np.zeros(nFull), units='m', desc='diameter along tower')

        # extra mass
        self.add_input('mass', 0.0, units='kg', desc='added mass')
        self.add_input('mI', np.zeros((6,)), units='kg*m**2', desc='mass moment of inertia about some point p [xx yy zz xy xz yz]')
        self.add_input('mrho', np.zeros((3,)), units='m', desc='xyz-location of p relative to node')
        self.add_input('transition_piece_mass', 0.0, units='kg', desc='point mass of transition piece')
        self.add_input('gravity_foundation_mass', 0.0, units='kg', desc='point mass of transition piece')
        self.add_input('transition_piece_height', 0.0, units='m', desc='height of transition piece above water line')
        self.add_input('foundation_height', 0.0, units='m', desc='height of foundation (0.0 for land, -water_depth for fixed bottom)')
        
        # point loads
        self.add_input('rna_F', np.zeros((3,)), units='N', desc='rna force')
        self.add_input('rna_M', np.zeros((3,)), units='N*m', desc='rna moment')

        # Monopile handling
        self.add_input('k_monopile', np.zeros(6), units='N/m', desc='Stiffness BCs for ocean soil.  Only used if monoflag inputis True')
        
        # spring reaction data.  Use float('inf') for rigid constraints.
        nK = nRefine+1 if self.options['monopile'] else 1
        self.add_output('kidx', np.zeros(nK, dtype=np.int_), desc='indices of z where external stiffness reactions should be applied.')
        self.add_output('kx', np.zeros(nK), units='N/m', desc='spring stiffness in x-direction')
        self.add_output('ky', np.zeros(nK), units='N/m', desc='spring stiffness in y-direction')
        self.add_output('kz', np.zeros(nK), units='N/m', desc='spring stiffness in z-direction')
        self.add_output('ktx', np.zeros(nK), units='N/m', desc='spring stiffness in theta_x-rotation')
        self.add_output('kty', np.zeros(nK), units='N/m', desc='spring stiffness in theta_y-rotation')
        self.add_output('ktz', np.zeros(nK), units='N/m', desc='spring stiffness in theta_z-rotation')
        
        # extra mass
        nMass = 3
        self.add_output('midx', np.zeros(nMass, dtype=np.int_), desc='indices where added mass should be applied.')
        self.add_output('m', np.zeros(nMass), units='kg', desc='added mass')
        self.add_output('mIxx', np.zeros(nMass), units='kg*m**2', desc='x mass moment of inertia about some point p')
        self.add_output('mIyy', np.zeros(nMass), units='kg*m**2', desc='y mass moment of inertia about some point p')
        self.add_output('mIzz', np.zeros(nMass), units='kg*m**2', desc='z mass moment of inertia about some point p')
        self.add_output('mIxy', np.zeros(nMass), units='kg*m**2', desc='xy mass moment of inertia about some point p')
        self.add_output('mIxz', np.zeros(nMass), units='kg*m**2', desc='xz mass moment of inertia about some point p')
        self.add_output('mIyz', np.zeros(nMass), units='kg*m**2', desc='yz mass moment of inertia about some point p')
        self.add_output('mrhox', np.zeros(nMass), units='m', desc='x-location of p relative to node')
        self.add_output('mrhoy', np.zeros(nMass), units='m', desc='y-location of p relative to node')
        self.add_output('mrhoz', np.zeros(nMass), units='m', desc='z-location of p relative to node')

        # point loads (if addGravityLoadForExtraMass=True be sure not to double count by adding those force here also)
        nPL = 1
        self.add_output('plidx', np.zeros(nPL, dtype=np.int_), desc='indices where point loads should be applied.')
        self.add_output('Fx', np.zeros(nPL), units='N', desc='point force in x-direction')
        self.add_output('Fy', np.zeros(nPL), units='N', desc='point force in y-direction')
        self.add_output('Fz', np.zeros(nPL), units='N', desc='point force in z-direction')
        self.add_output('Mxx', np.zeros(nPL), units='N*m', desc='point moment about x-axis')
        self.add_output('Myy', np.zeros(nPL), units='N*m', desc='point moment about y-axis')
        self.add_output('Mzz', np.zeros(nPL), units='N*m', desc='point moment about z-axis')

        #self.declare_partials('m','mass')
        #self.declare_partials(['mIxx','mIyy','mIzz','mIxy','mIxz','mIyz'], 'mI')
        #self.declare_partials(['Fx','Fy','Fz'], 'rna_F')
        #self.declare_partials(['Mxx','Myy','Mzz'], 'rna_M')

        
    def compute(self, inputs, outputs):
        nPoints = self.options['nPoints']
        nFull   = self.options['nFull']
        nRefine = int( (nFull-1)/(nPoints-1) )
        
        # Prepare RNA, transition piece, and gravity foundation (if any applicable) for "extra node mass"
        itrans = find_nearest(inputs['z'], inputs['transition_piece_height'])
        rtrans = 0.5*inputs['d'][itrans]
        mtrans = inputs['transition_piece_mass']
        Itrans = mtrans*rtrans**2. * np.r_[0.5, 0.5, 1.0, np.zeros(3)] # shell
        rgrav  = 0.5*inputs['d'][0]
        mgrav  = inputs['gravity_foundation_mass']
        Igrav  = mgrav*rgrav**2. * np.r_[0.25, 0.25, 0.5, np.zeros(3)] # disk
        # Note, need len()-1 because Frame3DD crashes if mass add at end
        outputs['midx']  = np.array([ len(inputs['z'])-1, itrans, 0 ], dtype=np.int_)
        outputs['m']     = np.array([ inputs['mass'],  mtrans,   mgrav ]).flatten()
        outputs['mIxx']  = np.array([ inputs['mI'][0], Itrans[0], Igrav[0] ]).flatten()
        outputs['mIyy']  = np.array([ inputs['mI'][1], Itrans[1], Igrav[1] ]).flatten()
        outputs['mIzz']  = np.array([ inputs['mI'][2], Itrans[2], Igrav[2] ]).flatten()
        outputs['mIxy']  = np.array([ inputs['mI'][3], Itrans[3], Igrav[3] ]).flatten()
        outputs['mIxz']  = np.array([ inputs['mI'][4], Itrans[4], Igrav[4] ]).flatten()
        outputs['mIyz']  = np.array([ inputs['mI'][5], Itrans[5], Igrav[5] ]).flatten()
        outputs['mrhox'] = np.array([ inputs['mrho'][0], 0.0, 0.0 ]).flatten()
        outputs['mrhoy'] = np.array([ inputs['mrho'][1], 0.0, 0.0 ]).flatten()
        outputs['mrhoz'] = np.array([ inputs['mrho'][2], 0.0, 0.0 ]).flatten()

        # Prepare point forces at RNA node
        outputs['plidx'] = np.array([ len(inputs['z'])-1 ], dtype=np.int_) # -1 b/c same reason as above
        outputs['Fx']    = np.array([ inputs['rna_F'][0] ]).flatten()
        outputs['Fy']    = np.array([ inputs['rna_F'][1] ]).flatten()
        outputs['Fz']    = np.array([ inputs['rna_F'][2] ]).flatten()
        outputs['Mxx']   = np.array([ inputs['rna_M'][0] ]).flatten()
        outputs['Myy']   = np.array([ inputs['rna_M'][1] ]).flatten()
        outputs['Mzz']   = np.array([ inputs['rna_M'][2] ]).flatten()

        # Prepare for reactions: rigid at tower base
        if self.options['monopile']:
            # Based on discretization, always same number of points in pile section
            kmono  = inputs['k_monopile']
            indvec = np.arange(0, nRefine+1, dtype=np.int_)
            outputs['kidx'] = indvec
            outputs['kx']   = kmono[0]*np.ones(indvec.shape)
            outputs['ky']   = kmono[2]*np.ones(indvec.shape)
            outputs['kz']   = kmono[4]*np.ones(indvec.shape)
            outputs['ktx']  = kmono[1]*np.ones(indvec.shape)
            outputs['kty']  = kmono[3]*np.ones(indvec.shape)
            outputs['ktz']  = kmono[5]*np.ones(indvec.shape)
        else:
            outputs['kidx'] = np.array([ 0 ], dtype=np.int_)
            outputs['kx']   = np.array([ 1.e16 ])
            outputs['ky']   = np.array([ 1.e16 ])
            outputs['kz']   = np.array([ 1.e16 ])
            outputs['ktx']  = np.array([ 1.e16 ])
            outputs['kty']  = np.array([ 1.e16 ])
            outputs['ktz']  = np.array([ 1.e16 ])

    '''
    def compute_partials(self, inputs, J, discrete_inputs):
        
        J['m','mass']    = 1.0
        J['mIxx','mI']   = np.eye(6)[0,:]
        J['mIyy','mI']   = np.eye(6)[1,:]
        J['mIzz','mI']   = np.eye(6)[2,:]
        J['mIxy','mI']   = np.eye(6)[3,:]
        J['mIxz','mI']   = np.eye(6)[4,:]
        J['mIyz','mI']   = np.eye(6)[5,:]
        J['Fx','rna_F']  = np.eye(3)[0,:]
        J['Fy','rna_F']  = np.eye(3)[2,:]
        J['Fz','rna_F']  = np.eye(3)[2,:]
        J['Mxx','rna_M'] = np.eye(3)[0,:]
        J['Myy','rna_M'] = np.eye(3)[2,:]
        J['Mzz','rna_M'] = np.eye(3)[2,:]
    '''

        
class TowerPostFrame(ExplicitComponent):
    def initialize(self):
        self.options.declare('nFull')
        #self.options.declare('nDEL')

    def setup(self):
        nFull = self.options['nFull']
        #nDEL  = self.options['nDEL']

        # effective geometry -- used for handbook methods to estimate hoop stress, buckling, fatigue
        self.add_input('z', np.zeros(nFull), units='m', desc='location along tower. start at bottom and go to top')
        self.add_input('d', np.zeros(nFull), units='m', desc='effective tower diameter for section')
        self.add_input('t', np.zeros(nFull-1), units='m', desc='effective shell thickness for section')
        self.add_input('L_reinforced', 0.0, units='m', desc='buckling length')

        # Material properties
        self.add_input('E', 0.0, units='N/m**2', desc='modulus of elasticity')

        # Processed Frame3DD outputs
        self.add_input('Fz', np.zeros(nFull-1), units='N', desc='Axial foce in vertical z-direction in cylinder structure.')
        self.add_input('Mxx', np.zeros(nFull-1), units='N*m', desc='Moment about x-axis in cylinder structure.')
        self.add_input('Myy', np.zeros(nFull-1), units='N*m', desc='Moment about y-axis in cylinder structure.')
        self.add_input('axial_stress', val=np.zeros(nFull-1), units='N/m**2', desc='axial stress in tower elements')
        self.add_input('shear_stress', val=np.zeros(nFull-1), units='N/m**2', desc='shear stress in tower elements')
        self.add_input('hoop_stress' , val=np.zeros(nFull-1), units='N/m**2', desc='hoop stress in tower elements')
        self.add_input('top_deflection_in', 0.0, units='m', desc='Deflection of tower top in yaw-aligned +x direction')

        # safety factors
        self.add_input('gamma_f', 1.35, desc='safety factor on loads')
        self.add_input('gamma_m', 1.1, desc='safety factor on materials')
        self.add_input('gamma_n', 1.0, desc='safety factor on consequence of failure')
        self.add_input('gamma_b', 1.1, desc='buckling safety factor')
        self.add_input('sigma_y', 0.0, units='N/m**2', desc='yield stress')
        self.add_input('gamma_fatigue', 1.755, desc='total safety factor for fatigue')

        # fatigue parameters
        self.add_input('life', 20.0, desc='fatigue life of tower')
        self.add_input('m_SN', 4, desc='slope of S/N curve')
        self.add_input('DC', 80.0, desc='standard value of stress')
        #self.add_input('z_DEL', np.zeros(nDEL), units='m', desc='absolute z coordinates of corresponding fatigue parameters')
        #self.add_input('M_DEL', np.zeros(nDEL), desc='fatigue parameters at corresponding z coordinates')

        # Frequencies
        self.add_input('f1', 0.0, units='Hz', desc='First natural frequency')
        self.add_input('f2', 0.0, units='Hz', desc='Second natural frequency')
        
        # outputs
        self.add_output('structural_frequencies', np.zeros(2), units='Hz', desc='First and second natural frequency')
        self.add_output('top_deflection', 0.0, units='m', desc='Deflection of tower top in yaw-aligned +x direction')
        self.add_output('stress', np.zeros(nFull-1), desc='Von Mises stress utilization along tower at specified locations.  incudes safety factor.')
        self.add_output('shell_buckling', np.zeros(nFull-1), desc='Shell buckling constraint.  Should be < 1 for feasibility.  Includes safety factors')
        self.add_output('global_buckling', np.zeros(nFull-1), desc='Global buckling constraint.  Should be < 1 for feasibility.  Includes safety factors')
        #self.add_output('damage', np.zeros(nFull-1), desc='Fatigue damage at each tower section')
        self.add_output('turbine_F', val=np.zeros(3), units='N', desc='Total force on tower+rna')
        self.add_output('turbine_M', val=np.zeros(3), units='N*m', desc='Total x-moment on tower+rna measured at base')
        
        # Derivatives
        # self.declare_partials('*', '*', method='fd', form='central', step=1e-6)

        
    def compute(self, inputs, outputs):
        # Unpack some variables
        axial_stress = inputs['axial_stress']
        shear_stress = inputs['shear_stress']
        hoop_stress  = inputs['hoop_stress']
        sigma_y      = inputs['sigma_y'] * np.ones(axial_stress.shape)
        E            = inputs['E'] * np.ones(axial_stress.shape)
        L_reinforced = inputs['L_reinforced'] * np.ones(axial_stress.shape)
        d,_          = nodal2sectional(inputs['d'])
        z_section,_  = nodal2sectional(inputs['z'])

        # Frequencies
        outputs['structural_frequencies'] = np.zeros(2)
        outputs['structural_frequencies'][0] = inputs['f1']
        outputs['structural_frequencies'][1] = inputs['f2']

        # Deflection
        outputs['top_deflection'] = inputs['top_deflection_in']
        
        # von mises stress
        outputs['stress'] = Util.vonMisesStressUtilization(axial_stress, hoop_stress, shear_stress,
                      inputs['gamma_f']*inputs['gamma_m']*inputs['gamma_n'], sigma_y)

        # shell buckling
        outputs['shell_buckling'] = Util.shellBucklingEurocode(d, inputs['t'], axial_stress, hoop_stress,
                                                                shear_stress, L_reinforced, E, sigma_y, inputs['gamma_f'], inputs['gamma_b'])

        # global buckling
        tower_height = inputs['z'][-1] - inputs['z'][0]
        M = np.sqrt(inputs['Mxx']**2 + inputs['Myy']**2)
        outputs['global_buckling'] = Util.bucklingGL(d, inputs['t'], inputs['Fz'], M, tower_height, E,
                                                      sigma_y, inputs['gamma_f'], inputs['gamma_b'])

        # fatigue
        N_DEL = 365.0*24.0*3600.0*inputs['life'] * np.ones(len(inputs['t']))
        #outputs['damage'] = np.zeros(N_DEL.shape)

        #if any(inputs['M_DEL']):
        #    M_DEL = np.interp(z_section, inputs['z_DEL'], inputs['M_DEL'])

        #    outputs['damage'] = Util.fatigue(M_DEL, N_DEL, d, inputs['t'], inputs['m_SN'],
        #                                      inputs['DC'], inputs['gamma_fatigue'], stress_factor=1.0, weld_factor=True)

# -----------------
#  Assembly
# -----------------

class TowerLeanSE(Group):

    def initialize(self):
        self.options.declare('nPoints')
        self.options.declare('nFull')
        self.options.declare('topLevelFlag', default=True)
        self.options.declare('monopile', default=False)
        
    def setup(self):
        nPoints      = self.options['nPoints']
        nFull        = self.options['nFull']
        topLevelFlag = self.options['topLevelFlag']
        nRefine      = int( (nFull-1)/(nPoints-1) )
        monopile     = self.options['monopile']
        
        # Independent variables that are unique to TowerSE
        towerIndeps = IndepVarComp()
        towerIndeps.add_output('tower_outer_diameter', np.zeros(nPoints), units='m')
        towerIndeps.add_output('tower_section_height', np.zeros(nPoints-1), units='m')
        towerIndeps.add_output('tower_wall_thickness', np.zeros(nPoints-1), units='m')
        towerIndeps.add_output('tower_buckling_length', 0.0, units='m')
        towerIndeps.add_output('tower_outfitting_factor', 0.0)
        towerIndeps.add_output('gravity_foundation_mass', 0.0, units='kg')
        towerIndeps.add_output('transition_piece_mass', 0.0, units='kg')
        towerIndeps.add_output('transition_piece_height', 0.0, units='m')
        towerIndeps.add_output('suctionpile_depth', 0.0, units='m')
        towerIndeps.add_output('suctionpile_depth_diam_ratio', 0.0)
        self.add_subsystem('towerIndepsLean', towerIndeps, promotes=['*'])

        # Independent variables that may be duplicated at higher levels of aggregation
        if topLevelFlag:
            sharedIndeps = IndepVarComp()
            sharedIndeps.add_output('hub_height', 0.0, units='m')
            sharedIndeps.add_output('material_density', 0.0, units='kg/m**3')
            self.add_subsystem('sharedIndepsLean', sharedIndeps, promotes=['*'])
        
        # All the static components
        self.add_subsystem('predisc', MonopileFoundation(monopile=monopile, nPoints=nPoints), promotes=['*'])
        # If doing fixed bottom monopile, we add an additional point for the pile (even for gravity foundations)
        if monopile:
            nPoints += 1
            nFull   += nRefine
            
        # Promote all but foundation_height so that we can override
        self.add_subsystem('geometry', CylinderDiscretization(nPoints=nPoints, nRefine=nRefine), promotes=['section_height','wall_thickness','diameter',
                                                                                                           'z_param','z_full','d_full','t_full'])
        
        self.add_subsystem('tgeometry', TowerDiscretization(nPoints=nPoints), promotes=['*'])
        
        self.add_subsystem('cm', CylinderMass(nPoints=nFull), promotes=['material_density','z_full','d_full','t_full',
                                                                        'material_cost_rate','labor_cost_rate','painting_cost_rate'])
        self.add_subsystem('tm', TowerMass(nFull=nFull), promotes=['z_full',
                                                                   'tower_mass','tower_center_of_mass','tower_section_center_of_mass','tower_I_base',
                                                                   'tower_raw_cost','gravity_foundation_mass','foundation_height',
                                                                   'transition_piece_mass','transition_piece_height',
                                                                   'monopile_mass','monopile_cost','monopile_length'])
        self.add_subsystem('gc', Util.GeometricConstraints(nPoints=nPoints), promotes=['min_d_to_t','max_taper','manufacturability','weldability','slope'])
        self.add_subsystem('turb', TurbineMass(), promotes=['turbine_mass','monopile_mass',
                                                            'tower_mass','tower_center_of_mass','tower_I_base',
                                                            'rna_mass', 'rna_cg', 'rna_I','hub_height'])
        
        # Connections for geometry and mass
        self.connect('foundation_height_out', 'geometry.foundation_height')
        self.connect('section_height_out', 'section_height')
        self.connect('outer_diameter_out', ['diameter', 'gc.d'])
        self.connect('wall_thickness_out', ['wall_thickness', 'gc.t'])
        self.connect('tower_outfitting_factor', 'cm.outfitting_factor')

        self.connect('cm.mass', 'tm.cylinder_mass')
        self.connect('cm.cost', 'tm.cylinder_cost')
        self.connect('cm.center_of_mass', 'tm.cylinder_center_of_mass')
        self.connect('cm.section_center_of_mass','tm.cylinder_section_center_of_mass')
        self.connect('cm.I_base','tm.cylinder_I_base')

        
class TowerSE(Group):

    def initialize(self):
        self.options.declare('nLC')
        self.options.declare('nPoints')
        self.options.declare('nFull')
        #self.options.declare('nDEL')
        self.options.declare('wind', default='')
        self.options.declare('topLevelFlag', default=True)
        self.options.declare('monopile', default=False)
    
    def setup(self):
        nLC           = self.options['nLC']
        nPoints       = self.options['nPoints']
        nFull         = self.options['nFull']
        # nDEL          = self.options['nDEL']
        wind          = self.options['wind']
        topLevelFlag  = self.options['topLevelFlag']
        monopile      = self.options['monopile']
        nRefine       = int( (nFull-1)/(nPoints-1) )
        nK = nRefine+1 if self.options['monopile'] else 1
        
        # Independent variables that are unique to TowerSE
        towerIndeps = IndepVarComp()
        towerIndeps.add_output('foundation_height', 0.0, units='m')
        #towerIndeps.add_output('tower_M_DEL', np.zeros(nDEL))
        #towerIndeps.add_output('tower_z_DEL', np.zeros(nDEL), units='m')
        towerIndeps.add_output('tower_force_discretization', 5.0)
        towerIndeps.add_output('soil_G', 0.0, units='N/m**2')
        towerIndeps.add_output('soil_nu', 0.0)
        towerIndeps.add_discrete_output('tower_add_gravity', True)
        towerIndeps.add_discrete_output('shear', True)
        towerIndeps.add_discrete_output('geom', True)
        towerIndeps.add_discrete_output('nM', 2)
        towerIndeps.add_discrete_output('Mmethod', 1)
        towerIndeps.add_discrete_output('lump', 0)
        towerIndeps.add_output('tol', 1e-9)
        towerIndeps.add_output('shift', 0.0)
        towerIndeps.add_output('DC', 0.0)
        towerIndeps.add_output('water_density', 1025.0, units='kg/m**3')
        towerIndeps.add_output('water_viscosity', 8.9e-4, units='kg/m/s')
        towerIndeps.add_output('significant_wave_height', 0.0, units='m')
        towerIndeps.add_output('significant_wave_period', 0.0, units='s')
        towerIndeps.add_output('wave_beta', 0.0, units='deg')
        self.add_subsystem('towerIndeps', towerIndeps, promotes=['*'])

        # Independent variables that may be duplicated at higher levels of aggregation
        if topLevelFlag:
            sharedIndeps = IndepVarComp()
            sharedIndeps.add_output('air_density', 1.225, units='kg/m**3')
            sharedIndeps.add_output('air_viscosity', 1.81206e-5, units='kg/m/s')
            sharedIndeps.add_output('shearExp', 0.0)
            sharedIndeps.add_output('wind_reference_height', 0.0, units='m')
            sharedIndeps.add_output('wind_z0', 0.0, units='m')
            sharedIndeps.add_output('wind_beta', 0.0, units='deg')
            sharedIndeps.add_output('cd_usr', -1.)
            sharedIndeps.add_output('yaw', 0.0, units='deg')
            sharedIndeps.add_output('E', 0.0, units='N/m**2')
            sharedIndeps.add_output('G', 0.0, units='N/m**2')
            sharedIndeps.add_output('sigma_y', 0.0, units='N/m**2')
            sharedIndeps.add_output('rna_mass', 0.0, units='kg')
            sharedIndeps.add_output('rna_cg', np.zeros(3), units='m')
            sharedIndeps.add_output('rna_I', np.zeros(6), units='kg*m**2')
            sharedIndeps.add_output('gamma_f', 0.0)
            sharedIndeps.add_output('gamma_m', 0.0)
            sharedIndeps.add_output('gamma_n', 0.0)
            sharedIndeps.add_output('gamma_b', 0.0)
            sharedIndeps.add_output('gamma_fatigue', 0.0)
            sharedIndeps.add_output('life', 0.0)
            sharedIndeps.add_output('m_SN', 0.0)
            self.add_subsystem('sharedIndeps', sharedIndeps, promotes=['*'])

        self.add_subsystem('geom', TowerLeanSE(nPoints=nPoints, nFull=nFull, topLevelFlag=topLevelFlag, monopile=monopile), promotes=['*'])
        # If doing fixed bottom monopile, we add an additional point for the pile (even for gravity foundations)
        if monopile:
            nPoints += 1
            nFull   += nRefine
        self.add_subsystem('props', CylindricalShellProperties(nFull=nFull))
        self.add_subsystem('soil', TowerSoil())

        # Connections for geometry and mass
        self.connect('d_full', 'props.d')
        self.connect('t_full', 'props.t')
        self.connect('d_full', 'soil.d0', src_indices=[0])
        self.connect('suctionpile_depth', 'soil.depth')
        self.connect('soil_G', 'soil.G')
        self.connect('soil_nu', 'soil.nu')
        
        # Add in all Components that drive load cases
        # Note multiple load cases have to be handled by replicating components and not groups/assemblies.
        # Replicating Groups replicates the IndepVarComps which doesn't play nicely in OpenMDAO
        for iLC in range(nLC):
            lc = '' if nLC==1 else str(iLC+1)
            
            if wind is None or wind.lower() in ['power', 'powerwind', '']:
                self.add_subsystem('wind'+lc, PowerWind(nPoints=nFull))
            elif wind.lower() == 'logwind':
                self.add_subsystem('wind'+lc, LogWind(nPoints=nFull))
            else:
                raise ValueError('Unknown wind type, '+wind)

            self.add_subsystem('windLoads'+lc, CylinderWindDrag(nPoints=nFull), promotes=['cd_usr'])

            if monopile:
                self.add_subsystem('wave'+lc, LinearWaves(nPoints=nFull), promotes=['z_floor'])
                self.add_subsystem('waveLoads'+lc, CylinderWaveDrag(nPoints=nFull), promotes=['cm','cd_usr'])

            self.add_subsystem('distLoads'+lc, AeroHydroLoads(nPoints=nFull), promotes=['yaw'])

            self.add_subsystem('pre'+lc, TowerPreFrame(nFull=nFull, nPoints=nPoints, monopile=monopile), promotes=['transition_piece_mass',
                                                                                                                   'transition_piece_height',
                                                                                                                   'gravity_foundation_mass'])
            self.add_subsystem('tower'+lc, CylinderFrame3DD(npts=nFull, nK=nK, nMass=3, nPL=1), promotes=['E','G','tol','Mmethod','geom','lump','shear',
                                                                                                         'nM','shift','sigma_y'])
            self.add_subsystem('post'+lc, TowerPostFrame(nFull=nFull), promotes=['E','sigma_y','DC','life','m_SN',
                                                                                 'gamma_b','gamma_f','gamma_fatigue','gamma_m','gamma_n'])
            
            self.connect('z_full', ['wind'+lc+'.z', 'windLoads'+lc+'.z', 'distLoads'+lc+'.z', 'pre'+lc+'.z', 'tower'+lc+'.z', 'post'+lc+'.z'])
            self.connect('d_full', ['windLoads'+lc+'.d', 'tower'+lc+'.d', 'pre'+lc+'.d', 'post'+lc+'.d'])
            if monopile:
                self.connect('z_full', ['wave'+lc+'.z', 'waveLoads'+lc+'.z'])
                self.connect('d_full', 'waveLoads'+lc+'.d')

            if topLevelFlag:
                self.connect('rna_mass', 'pre'+lc+'.mass')
                self.connect('rna_cg', 'pre'+lc+'.mrho')
                self.connect('rna_I', 'pre'+lc+'.mI')
                self.connect('material_density', 'tower'+lc+'.rho')

            self.connect('pre'+lc+'.kidx', 'tower'+lc+'.kidx')
            self.connect('pre'+lc+'.kx', 'tower'+lc+'.kx')
            self.connect('pre'+lc+'.ky', 'tower'+lc+'.ky')
            self.connect('pre'+lc+'.kz', 'tower'+lc+'.kz')
            self.connect('pre'+lc+'.ktx', 'tower'+lc+'.ktx')
            self.connect('pre'+lc+'.kty', 'tower'+lc+'.kty')
            self.connect('pre'+lc+'.ktz', 'tower'+lc+'.ktz')
            self.connect('pre'+lc+'.midx', 'tower'+lc+'.midx')
            self.connect('pre'+lc+'.m', 'tower'+lc+'.m')
            self.connect('pre'+lc+'.mIxx', 'tower'+lc+'.mIxx')
            self.connect('pre'+lc+'.mIyy', 'tower'+lc+'.mIyy')
            self.connect('pre'+lc+'.mIzz', 'tower'+lc+'.mIzz')
            self.connect('pre'+lc+'.mIxy', 'tower'+lc+'.mIxy')
            self.connect('pre'+lc+'.mIxz', 'tower'+lc+'.mIxz')
            self.connect('pre'+lc+'.mIyz', 'tower'+lc+'.mIyz')
            self.connect('pre'+lc+'.mrhox', 'tower'+lc+'.mrhox')
            self.connect('pre'+lc+'.mrhoy', 'tower'+lc+'.mrhoy')
            self.connect('pre'+lc+'.mrhoz', 'tower'+lc+'.mrhoz')

            self.connect('pre'+lc+'.plidx', 'tower'+lc+'.plidx')
            self.connect('pre'+lc+'.Fx', 'tower'+lc+'.Fx')
            self.connect('pre'+lc+'.Fy', 'tower'+lc+'.Fy')
            self.connect('pre'+lc+'.Fz', 'tower'+lc+'.Fz')
            self.connect('pre'+lc+'.Mxx', 'tower'+lc+'.Mxx')
            self.connect('pre'+lc+'.Myy', 'tower'+lc+'.Myy')
            self.connect('pre'+lc+'.Mzz', 'tower'+lc+'.Mzz')
            self.connect('tower_force_discretization', 'tower'+lc+'.dx')
            self.connect('tower_add_gravity', 'tower'+lc+'.addGravityLoadForExtraMass')
            self.connect('t_full', ['tower'+lc+'.t','post'+lc+'.t'])
            self.connect('soil.k', 'pre'+lc+'.k_monopile')

            self.connect('tower'+lc+'.f1', 'post'+lc+'.f1')
            self.connect('tower'+lc+'.f2', 'post'+lc+'.f2')
            self.connect('tower'+lc+'.Fz_out', 'post'+lc+'.Fz')
            self.connect('tower'+lc+'.Mxx_out', 'post'+lc+'.Mxx')
            self.connect('tower'+lc+'.Myy_out', 'post'+lc+'.Myy')
            self.connect('tower'+lc+'.axial_stress', 'post'+lc+'.axial_stress')
            self.connect('tower'+lc+'.shear_stress', 'post'+lc+'.shear_stress')
            self.connect('tower'+lc+'.hoop_stress_euro', 'post'+lc+'.hoop_stress')
            self.connect('tower'+lc+'.top_deflection', 'post'+lc+'.top_deflection_in')
        
            # connections to wind, wave
            self.connect('wind'+lc+'.U', 'windLoads'+lc+'.U')
            if topLevelFlag:
                self.connect('wind_reference_height', 'wind'+lc+'.zref')
                self.connect('wind_z0', 'wind'+lc+'.z0')
                if monopile:
                    self.connect('wind_z0', 'wave'+lc+'.z_surface')
                #self.connect('z_floor', 'waveLoads'+lc+'.wlevel')
                if wind=='PowerWind':
                    self.connect('shearExp', 'wind'+lc+'.shearExp')
                
                # connections to windLoads1
                self.connect('air_density', 'windLoads'+lc+'.rho')
                self.connect('air_viscosity', 'windLoads'+lc+'.mu')
                self.connect('wind_beta', 'windLoads'+lc+'.beta')

            if monopile:
                # connections to waveLoads1
                self.connect('water_density', ['wave'+lc+'.rho', 'waveLoads'+lc+'.rho'])
                self.connect('water_viscosity', 'waveLoads'+lc+'.mu')
                self.connect('wave_beta', 'waveLoads'+lc+'.beta')
                self.connect('significant_wave_height', 'wave'+lc+'.hmax')
                self.connect('significant_wave_period', 'wave'+lc+'.T')
                self.connect('foundation_height', 'z_floor')
                    
                self.connect('wave'+lc+'.U', 'waveLoads'+lc+'.U')
                self.connect('wave'+lc+'.A', 'waveLoads'+lc+'.A')
                self.connect('wave'+lc+'.p', 'waveLoads'+lc+'.p')

            # connections to distLoads1
            self.connect('windLoads'+lc+'.windLoads_Px', 'distLoads'+lc+'.windLoads_Px')
            self.connect('windLoads'+lc+'.windLoads_Py', 'distLoads'+lc+'.windLoads_Py')
            self.connect('windLoads'+lc+'.windLoads_Pz', 'distLoads'+lc+'.windLoads_Pz')
            self.connect('windLoads'+lc+'.windLoads_qdyn', 'distLoads'+lc+'.windLoads_qdyn')
            self.connect('windLoads'+lc+'.windLoads_beta', 'distLoads'+lc+'.windLoads_beta')
            #self.connect('windLoads'+lc+'.windLoads_Px0', 'distLoads'+lc+'.windLoads_Px0')
            #self.connect('windLoads'+lc+'.windLoads_Py0', 'distLoads'+lc+'.windLoads_Py0')
            #self.connect('windLoads'+lc+'.windLoads_Pz0', 'distLoads'+lc+'.windLoads_Pz0')
            #self.connect('windLoads'+lc+'.windLoads_qdyn0', 'distLoads'+lc+'.windLoads_qdyn0')
            #self.connect('windLoads'+lc+'.windLoads_beta0', 'distLoads'+lc+'.windLoads_beta0')
            self.connect('windLoads'+lc+'.windLoads_z', 'distLoads'+lc+'.windLoads_z')
            self.connect('windLoads'+lc+'.windLoads_d', 'distLoads'+lc+'.windLoads_d')

            if monopile:
                self.connect('waveLoads'+lc+'.waveLoads_Px', 'distLoads'+lc+'.waveLoads_Px')
                self.connect('waveLoads'+lc+'.waveLoads_Py', 'distLoads'+lc+'.waveLoads_Py')
                self.connect('waveLoads'+lc+'.waveLoads_Pz', 'distLoads'+lc+'.waveLoads_Pz')
                self.connect('waveLoads'+lc+'.waveLoads_pt', 'distLoads'+lc+'.waveLoads_qdyn')
                self.connect('waveLoads'+lc+'.waveLoads_beta', 'distLoads'+lc+'.waveLoads_beta')
                #self.connect('waveLoads'+lc+'.waveLoads_Px0', 'distLoads'+lc+'.waveLoads_Px0')
                #self.connect('waveLoads'+lc+'.waveLoads_Py0', 'distLoads'+lc+'.waveLoads_Py0')
                #self.connect('waveLoads'+lc+'.waveLoads_Pz0', 'distLoads'+lc+'.waveLoads_Pz0')
                #self.connect('waveLoads'+lc+'.waveLoads_qdyn0', 'distLoads'+lc+'.waveLoads_qdyn0')
                #self.connect('waveLoads'+lc+'.waveLoads_beta0', 'distLoads'+lc+'.waveLoads_beta0')
                self.connect('waveLoads'+lc+'.waveLoads_z', 'distLoads'+lc+'.waveLoads_z')
                self.connect('waveLoads'+lc+'.waveLoads_d', 'distLoads'+lc+'.waveLoads_d')

            # Tower connections
            self.connect('tower_buckling_length', ['tower'+lc+'.L_reinforced', 'post'+lc+'.L_reinforced'])
            #self.connect('tower_M_DEL', 'post'+lc+'.M_DEL')
            #self.connect('tower_z_DEL', 'post'+lc+'.z_DEL')

            self.connect('props.Az', 'tower'+lc+'.Az')
            self.connect('props.Asx', 'tower'+lc+'.Asx')
            self.connect('props.Asy', 'tower'+lc+'.Asy')
            self.connect('props.Jz', 'tower'+lc+'.Jz')
            self.connect('props.Ixx', 'tower'+lc+'.Ixx')
            self.connect('props.Iyy', 'tower'+lc+'.Iyy')

            self.connect('distLoads'+lc+'.Px',   'tower'+lc+'.Px')
            self.connect('distLoads'+lc+'.Py',   'tower'+lc+'.Py')
            self.connect('distLoads'+lc+'.Pz',   'tower'+lc+'.Pz')
            self.connect('distLoads'+lc+'.qdyn', 'tower'+lc+'.qdyn')

        
if __name__ == '__main__':
    # --- tower setup ------
    from wisdem.commonse.environment import PowerWind
    from wisdem.commonse.environment import LogWind

    optimize_flag = False

    # --- geometry ----
    h_param = np.diff(np.array([0.0, 43.8, 87.6]))
    d_param = np.array([6.0, 4.935, 3.87])
    t_param = 1.3*np.array([0.025, 0.021])
    z_foundation = 0.0
    L_reinforced = 30.0  # [m] buckling length
    theta_stress = 0.0
    yaw = 0.0
    Koutfitting = 1.07

    # --- material props ---
    E = 210e9
    G = 80.8e9
    rho = 8500.0
    sigma_y = 450.0e6

    # --- extra mass ----
    m = np.array([285598.8])
    mIxx = 1.14930678e+08
    mIyy = 2.20354030e+07
    mIzz = 1.87597425e+07
    mIxy = 0.0
    mIxz = 5.03710467e+05
    mIyz = 0.0
    mI = np.array([mIxx, mIyy, mIzz, mIxy, mIxz, mIyz])
    mrho = np.array([-1.13197635, 0.0, 0.50875268])
    # -----------

    # --- wind ---
    wind_zref = 90.0
    wind_z0 = 0.0
    shearExp = 0.2
    cd_usr = -1.
    # ---------------

    # --- wave ---
    hmax = 0.0
    T = 1.0
    cm = 1.0
    monopile = False
    suction_depth = 0.0
    soilG = 140e6
    soilnu = 0.4
    # ---------------

    
    # two load cases.  TODO: use a case iterator
    
    # # --- loading case 1: max Thrust ---
    wind_Uref1 = 11.73732
    Fx1 = 1284744.19620519
    Fy1 = 0.
    Fz1 = -2914124.84400512 + m*gravity
    Mxx1 = 3963732.76208099
    Myy1 = -2275104.79420872
    Mzz1 = -346781.68192839
    # # ---------------

    # # --- loading case 2: max wind speed ---
    wind_Uref2 = 70.0
    Fx2 = 930198.60063279
    Fy2 = 0.
    Fz2 = -2883106.12368949 + m*gravity
    Mxx2 = -1683669.22411597
    Myy2 = -2522475.34625363
    Mzz2 = 147301.97023764
    # # ---------------

    # --- safety factors ---
    gamma_f = 1.35
    gamma_m = 1.3
    gamma_n = 1.0
    gamma_b = 1.1
    # ---------------

    # --- fatigue ---
    z_DEL = np.array([0.000, 1.327, 3.982, 6.636, 9.291, 11.945, 14.600, 17.255, 19.909, 22.564, 25.218, 27.873, 30.527, 33.182, 35.836, 38.491, 41.145, 43.800, 46.455, 49.109, 51.764, 54.418, 57.073, 59.727, 62.382, 65.036, 67.691, 70.345, 73.000, 75.655, 78.309, 80.964, 83.618, 86.273, 87.600])
    M_DEL = 1e3*np.array([8.2940E+003, 8.1518E+003, 7.8831E+003, 7.6099E+003, 7.3359E+003, 7.0577E+003, 6.7821E+003, 6.5119E+003, 6.2391E+003, 5.9707E+003, 5.7070E+003, 5.4500E+003, 5.2015E+003, 4.9588E+003, 4.7202E+003, 4.4884E+003, 4.2577E+003, 4.0246E+003, 3.7942E+003, 3.5664E+003, 3.3406E+003, 3.1184E+003, 2.8977E+003, 2.6811E+003, 2.4719E+003, 2.2663E+003, 2.0673E+003, 1.8769E+003, 1.7017E+003, 1.5479E+003, 1.4207E+003, 1.3304E+003, 1.2780E+003, 1.2673E+003, 1.2761E+003])
    nDEL = len(z_DEL)
    gamma_fatigue = 1.35*1.3*1.0
    life = 20.0
    m_SN = 4
    # ---------------


    # --- constraints ---
    min_d_to_t   = 120.0
    max_taper    = 0.2
    # ---------------

    # # V_max = 80.0  # tip speed
    # # D = 126.0
    # # .freq1p = V_max / (D/2) / (2*pi)  # convert to Hz

    nPoints = len(d_param)
    nFull   = 3*(nPoints-1) + 1
    wind = 'PowerWind'
    nLC = 2
    
    prob = Problem()
    prob.model = TowerSE(nLC=nLC, nPoints=nPoints, nFull=nFull, wind=wind, topLevelFlag=True, monopile=monopile)

    if optimize_flag:

        # --- optimizer imports ---
        from openmdao.api import ScipyOptimizeDriver, SqliteRecorder
        # ----------------------

        # --- Setup Pptimizer ---
        prob.driver = ScipyOptimizeDriver()
        prob.driver.options['optimizer'] = 'SLSQP'
        # ----------------------

        # --- Objective ---
        prob.model.add_objective('tower_mass', scaler=1e-6)
        # ----------------------

        # --- Design Variables ---
        prob.model.add_design_var('tower_outer_diameter', lower=3.87, upper=6.0, indices=[m for m in range(nPoints-1)])
        prob.model.add_design_var('tower_wall_thickness', lower=4e-3, upper=2e-1)
        prob.model.add_design_var('tower_section_height', lower=2.0, upper=80.0)
        # ----------------------

        # --- recorder ---
        frecord_opt = 'tower_example-optimizer.sql'
        frecord = 'tower_example.sql'
        prob.driver.add_recorder(SqliteRecorder(frecord_opt))
        prob.add_recorder(SqliteRecorder(frecord))
        # ----------------------

        # --- Constraints ---
        prob.model.add_constraint('height_constraint',    lower=-1e-2,upper=1.e-2)
        prob.model.add_constraint('post1.stress',          upper=1.0)
        prob.model.add_constraint('post1.global_buckling', upper=1.0)
        prob.model.add_constraint('post1.shell_buckling',  upper=1.0)
        prob.model.add_constraint('post2.stress',          upper=1.0)
        prob.model.add_constraint('post2.global_buckling', upper=1.0)
        prob.model.add_constraint('post2.shell_buckling',  upper=1.0)
        prob.model.add_constraint('weldability',          upper=0.0)
        prob.model.add_constraint('manufacturability',    lower=0.0)
        prob.model.add_constraint('slope',                upper=1.0)
        prob.model.add_constraint('tower1.f1',             lower=0.3)
        # ----------------------

    prob.setup()

    if wind=='PowerWind':
        prob['shearExp'] = shearExp

    # assign values to params

    # --- geometry ----
    prob['hub_height'] = h_param.sum()
    prob['foundation_height'] = 0.0
    prob['tower_section_height'] = h_param
    prob['tower_outer_diameter'] = d_param
    prob['tower_wall_thickness'] = t_param
    prob['tower_buckling_length'] = L_reinforced
    prob['tower_outfitting_factor'] = Koutfitting
    prob['yaw'] = yaw
    prob['suctionpile_depth'] = suction_depth
    prob['suctionpile_depth_diam_ratio'] = 3.35
    prob['soil_G'] = soilG
    prob['soil_nu'] = soilnu
    # --- material props ---
    prob['E'] = E
    prob['G'] = G
    prob['material_density'] = rho
    prob['sigma_y'] = sigma_y

    # --- extra mass ----
    prob['rna_mass'] = m
    prob['rna_I'] = mI
    prob['rna_cg'] = mrho
    # -----------

    # --- wind & wave ---
    prob['wind_reference_height'] = wind_zref
    prob['wind_z0'] = wind_z0
    prob['cd_usr'] = cd_usr
    prob['air_density'] = 1.225
    prob['air_viscosity'] = 1.7934e-5
    prob['water_density'] = 1025.0
    prob['water_viscosity'] = 1.3351e-3
    prob['wind_beta'] = prob['wave_beta'] = 0.0
    prob['significant_wave_height'] = hmax
    prob['significant_wave_period'] = T
    #prob['waveLoads1.U0'] = prob['waveLoads1.A0'] = prob['waveLoads1.beta0'] = prob['waveLoads2.U0'] = prob['waveLoads2.A0'] = prob['waveLoads2.beta0'] = 0.0
    # ---------------

    # --- safety factors ---
    prob['gamma_f'] = gamma_f
    prob['gamma_m'] = gamma_m
    prob['gamma_n'] = gamma_n
    prob['gamma_b'] = gamma_b
    prob['gamma_fatigue'] = gamma_fatigue
    # ---------------

    prob['DC'] = 80.0
    prob['shear'] = True
    prob['geom'] = True
    prob['tower_force_discretization'] = 5.0
    prob['nM'] = 2
    prob['Mmethod'] = 1
    prob['lump'] = 0
    prob['tol'] = 1e-9
    prob['shift'] = 0.0

    
    # --- fatigue ---
    #prob['tower_z_DEL'] = z_DEL
    #prob['tower_M_DEL'] = M_DEL
    prob['life'] = life
    prob['m_SN'] = m_SN
    # ---------------

    # --- constraints ---
    prob['min_d_to_t'] = min_d_to_t
    prob['max_taper'] = max_taper
    # ---------------


    # # --- loading case 1: max Thrust ---
    prob['wind1.Uref'] = wind_Uref1

    prob['pre1.rna_F'] = np.array([Fx1, Fy1, Fz1])
    prob['pre1.rna_M'] = np.array([Mxx1, Myy1, Mzz1])
    # # ---------------


    # # --- loading case 2: max Wind Speed ---
    prob['wind2.Uref'] = wind_Uref2

    prob['pre2.rna_F'] = np.array([Fx2, Fy2, Fz2])
    prob['pre2.rna_M' ] = np.array([Mxx2, Myy2, Mzz2])

    # # --- run ---
    if optimize_flag:
        prob.model.approx_totals()
        prob.run_driver()
    else:
        prob.run_model()

    z,_ = nodal2sectional(prob['z_full'])

    print('zs=', prob['z_full'])
    print('ds=', prob['d_full'])
    print('ts=', prob['t_full'])
    print('mass (kg) =', prob['tower_mass'])
    print('cg (m) =', prob['tower_center_of_mass'])
    print('weldability =', prob['weldability'])
    print('manufacturability =', prob['manufacturability'])
    print('\nwind: ', prob['wind1.Uref'])
    print('f1 (Hz) =', prob['tower1.f1'])
    print('top_deflection1 (m) =', prob['post1.top_deflection'])
    print('stress1 =', prob['post1.stress'])
    print('GL buckling =', prob['post1.global_buckling'])
    print('Shell buckling =', prob['post1.shell_buckling'])
    #print('damage =', prob['post1.damage'])
    print('\nwind: ', prob['wind2.Uref'])
    print('f1 (Hz) =', prob['tower2.f1'])
    print('top_deflection2 (m) =', prob['post2.top_deflection'])
    print('stress2 =', prob['post2.stress'])
    print('GL buckling =', prob['post2.global_buckling'])
    print('Shell buckling =', prob['post2.shell_buckling'])
    #print('damage =', prob['post2.damage'])


    stress1 = np.copy( prob['post1.stress'] )
    shellBuckle1 = np.copy( prob['post1.shell_buckling'] )
    globalBuckle1 = np.copy( prob['post1.global_buckling'] )
    #damage1 = np.copy( prob['post1.damage'] )

    stress2 = prob['post2.stress']
    shellBuckle2 = prob['post2.shell_buckling']
    globalBuckle2 = prob['post2.global_buckling']
    #damage2 = prob['post2.damage']

    
    import matplotlib.pyplot as plt
    plt.figure(figsize=(5.0, 3.5))
    plt.subplot2grid((3, 3), (0, 0), colspan=2, rowspan=3)
    plt.plot(stress1, z, label='stress 1')
    plt.plot(stress2, z, label='stress 2')
    plt.plot(shellBuckle1, z, label='shell buckling 1')
    plt.plot(shellBuckle2, z, label='shell buckling 2')
    plt.plot(globalBuckle1, z, label='global buckling 1')
    plt.plot(globalBuckle2, z, label='global buckling 2')
    #plt.plot(damage1, z, label='damage 1')
    #plt.plot(damage2, z, label='damage 2')
    plt.legend(bbox_to_anchor=(1.05, 1.0), loc=2)
    plt.xlabel('utilization')
    plt.ylabel('height along tower (m)')

    #plt.figure(2)
    #plt.plot(prob['d_full']/2.+max(prob['d_full']), z, 'ok')
    #plt.plot(prob['d_full']/-2.+max(prob['d_full']), z, 'ok')

    #fig = plt.figure(3)
    #ax1 = fig.add_subplot(121)
    #ax2 = fig.add_subplot(122)

    #ax1.plot(prob['wind1.U'], z)
    #ax2.plot(prob['wind2.U'], z)
    #plt.tight_layout()
    plt.show()

    print(prob['tower1.base_F'])
    print(prob['tower1.base_M'])
    print(prob['tower2.base_F'])
    print(prob['tower2.base_M'])
    # ------------




    
    
