#(c) Rory Clune 05/16/2011

from numpy import*
from Structures import *
from SearchAlgorithms import*
import matplotlib.pyplot as plt

if __name__ == '__main__':
    #-----------
    # INPUT DATA FOR A REALLY SIMPLE BEAM STRUCTURE
    #-----------
    # The geometry of the structural nodes (10-bar truss)
    geometry = array([[0., 0. ], [ 5., 2. ], [ 10., 0. ]])
    # A connectivity matrix defines how nodes are connected with beams. a row contains ( beam - node - node )
    connectivity = array([[ 0, 1 ], [ 1, 2 ]])
    # Restrained degrees of freedom - pins at either end
    fixity = array([0,1,6,7])
    # The loading - a downward point load at midspan
    load = zeros((3*geometry.shape[0])) 
    load[4] = -100
    A = 10*ones(connectivity.shape[0])
    I = 10*ones(connectivity.shape[0])    
    s = structure(geometry, connectivity, fixity, load, A, I)
    print 'maximum displacement at start: ' + str(amax(fabs(s.analyze())))
    s.show(40,40,100,50,array([4]));    
    
    #-----------
    # RUN OPTIMIZATION
    #-----------    
    s.Objective = s.VolumeObjective
    opt_geom = geometry[1,:]
    opt_geom = append(opt_geom, (s.A, s.I))
    
    def input_x(self, x):
        self.geometry[1,:] = x[0:2]
        self.A = x[2:4]
        self.I = x[4:6]
    s.input_x = input_x
    
    UB = array([7, 10, 2000, 2000, 2000, 2000], dtype = float)
    LB = array([3, -10, 1, 1, 1, 1], dtype = float)
#    UB = array([7, 10], dtype = float)
#    LB = array([3, -10], dtype = float)
        
    scale = array([1, 1, 1, 1, 1, 1], dtype = float)
#    scale = array([1, 1], dtype = float)
    
    opt_geom *= (1/scale); UB *= (1/scale); LB *= (1/scale)    
    s.displacement_limit = 0.4
    
    optimizer = NelderMeadConstrained(opt_geom.shape[0], s, opt_geom, UB, LB, 1e-10, 1000)
    optimizer.s = 100 #controls the rate at which a violation of a constraint causes the multiplier to increase
    optimizer.scale = scale
    optimizer.penalty = 0*ones(3*geometry.shape[0]) 
    optimizer.Run()
    
    print 'maximum displacement at end: ' + str(amax(fabs(s.analyze())))
#    plt.plot(optimizer.obj_history); plt.title('Objective function value achieved'); plt.show()   
    
    s.show(40,40,100,50,array([4]));