# Path: aeropy/structural/stable_solution.py
# class properties():
#     def __init__(self, young=70e9, poisson=.3, dimensions=[0.01, 0.01],
#                  crosssection='square'):
#         """Stores calculates material or structural properties. If 'square',
#         dimensions = [width, height]"""
#         self.young = young
#         self.poisson = poisson
#         self.dimensions = dimensions
#         self.length = 1
# 
#         if crosssection == 'square':
#             self.area = dimensions[0]*dimensions[1]
#             self.inertia = (dimensions[0]*dimensions[1]**3)/12.
#         if young < 0 or poisson < 0 or poisson > 1 or self.area <= 0:
#             raise Exception('Material properties need to make sense')
# 
# class boundary_conditions():
#     def __init__(self, concentrated_load=np.array([[0, 0, 0], ]), load_x=[1],
#                        distributed_load = 0):
#         self.concentrated_load = concentrated_load
#         self.concentrated_x = load_x
#         self.concentrated_n = len(concentrated_load)
#         self.distributed_load = None
# 
#         self.distributed_load = distributed_load
# 
#         if len(concentrated_load) != len(load_x):
#             raise Exception('load values and x lists have to match')
# 
# class euler_bernoulle():
#     def __init__(self, properties, load, load_type, geometry,
#                  x = np.linspace(0,1), boundaries = {'x':[], 'u':[]}):
#         self.properties = properties
#         self.load = load
#         self.load_type = load_type
#         self.g = geometry
#         self.g.x1_grid = x
#         self.boundaries = boundaries
# 
#     def analytical_solutions(self):
#         if self.load_type == 'concentrated':
#             bp = self.properties
#             self.g.D = self.load/(6*bp.young*bp.inertia) * \
#                             np.array([0, 0, 3, -1, 0])
#         elif self.load_type == 'distributed':
#             bp = self.properties
#             c2 = 6*(bp.length**2)*self.load/(24*bp.young*bp.inertia)
#             c3 = -4*(bp.length)*self.load/(24*bp.young*bp.inertia)
#             c4 = (1)*self.load/(24*bp.young*bp.inertia)
#             self.g.D = np.array([0, 0, c2, c3, c4])
# 
#     def bending_strain(self):
#         print(self.g.x1_grid)
#         self.B = self.g.x3(self.g.x1_grid, 'x11')
# 
#     def free_energy(self):
#         bp = self.properties
#         ddu = self.g.x3(self.g.x1_grid, 'x11')
#         self.phi = bp.young*bp.inertia/2*ddu**2
# 
#     def strain_energy(self):
#         self.U = np.trapz(self.phi, self.g.x1_grid)
# 
#     def work(self):
#         u = self.g.x3(self.g.x1_grid)
#         if self.load_type == 'concentrated':
#             self.W = self.load*u[-1]
#         elif self.load_type == 'distributed':
#             self.W = np.trapz(self.load*u, self.g.x1_grid)
# 
#     def residual(self):
#         self.R = self.U - self.W
# 
#     def minimum_potential(self, x0=[0,0,0]):
#         def to_optimize(x):
#             self.g.D = [0,0] + list(x)
#             self.work()
#             self.free_energy()
#             self.strain_energy()
#             self.residual()
#             return self.R
# 
#         # With bounds
#         bounds = np.array(((-0.01,0.01),)*len(x0))
# 
#         res = minimize(to_optimize, x0, bounds=bounds)
#         self.g.D = [0,0] + list(res.x)
#         self.work()
#         self.free_energy()
#         self.strain_energy()
#         self.residual()
#         return(res.x, res.fun)


from aeropy.geometry.parametric import CoordinateSystem
from aeropy.structural.stable_solution import (properties, boundary_conditions, euler_bernoulle)
from aeropy.structural.shell import shell
from aeropy.xfoil_module import output_reader
from scipy.integrate import quad
from scipy.optimize import minimize
from optimization_tools.DOE import DOE
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np
import pickle


# beam.minimum_potential(x0=[0,0], input_function = lambda x: [0,0] + list(x),
                       # bounds = bounds)

# Euler Bernoulle
EB = euler_bernoulle(bp, -1, 'concentrated', beam.g_c, x=beam.g_c.x1_grid)
EB.analytical_solutions()
EB.bending_strain()
EB.free_energy()
EB.strain_energy()
print('U_EB: ', EB.U)
print('U_PL: ', beam.U)
# mod = 1/(1-0.3**2)
# print('U_corrected: ', beam.U/mod)
#
print('Bending_strain_EB: ', EB.B)
print('Bending_strain_PL: ', beam.rho)
#
# print('Phi_EB: ', EB.phi)
# print('Phi_PL: ', beam.phi)
# print('Phi_PL_B: ', beam.phi_B)
#
# print('ddx3_EB: ', EB.g.x3(beam.g_c.x1_grid, 'x11'))
# print('ddx3_PL: ', beam.g_c.x3(beam.g_c.x1_grid, 'x11'))

print('theta1',beam.g_c.x1(beam.g_c.x1_grid, 'theta1'))
print('theta11', beam.g_c.x1(beam.g_c.x1_grid, 'theta11'))
P = -1
E = bp.young

# Path: aeropy/structural/stable_solution.py
# class properties():
#     def __init__(self, young=70e9, poisson=.3, dimensions=[0.01, 0.01],
#                  crosssection='square'):
#         """Stores calculates material or structural properties. If 'square',
#         dimensions = [width, height]"""
#         self.young = young
#         self.poisson = poisson
#         self.dimensions = dimensions
#         self.length = 1
# 
#         if crosssection == 'square':
#             self.area = dimensions[0]*dimensions[1]
#             self.inertia = (dimensions[0]*dimensions[1]**3)/12.
#         if young < 0 or poisson < 0 or poisson > 1 or self.area <= 0:
#             raise Exception('Material properties need to make sense')
# 
# class boundary_conditions():
#     def __init__(self, concentrated_load=np.array([[0, 0, 0], ]), load_x=[1],
#                        distributed_load = 0):
#         self.concentrated_load = concentrated_load
#         self.concentrated_x = load_x
#         self.concentrated_n = len(concentrated_load)
#         self.distributed_load = None
# 
#         self.distributed_load = distributed_load
# 
#         if len(concentrated_load) != len(load_x):
#             raise Exception('load values and x lists have to match')
# 
# class euler_bernoulle():
#     def __init__(self, properties, load, load_type, geometry,
#                  x = np.linspace(0,1), boundaries = {'x':[], 'u':[]}):
#         self.properties = properties
#         self.load = load
#         self.load_type = load_type
#         self.g = geometry
#         self.g.x1_grid = x
#         self.boundaries = boundaries
# 
#     def analytical_solutions(self):
#         if self.load_type == 'concentrated':
#             bp = self.properties
#             self.g.D = self.load/(6*bp.young*bp.inertia) * \
#                             np.array([0, 0, 3, -1, 0])
#         elif self.load_type == 'distributed':
#             bp = self.properties
#             c2 = 6*(bp.length**2)*self.load/(24*bp.young*bp.inertia)
#             c3 = -4*(bp.length)*self.load/(24*bp.young*bp.inertia)
#             c4 = (1)*self.load/(24*bp.young*bp.inertia)
#             self.g.D = np.array([0, 0, c2, c3, c4])
# 
#     def bending_strain(self):
#         print(self.g.x1_grid)
#         self.B = self.g.x3(self.g.x1_grid, 'x11')
# 
#     def free_energy(self):
#         bp = self.properties
#         ddu = self.g.x3(self.g.x1_grid, 'x11')
#         self.phi = bp.young*bp.inertia/2*ddu**2
# 
#     def strain_energy(self):
#         self.U = np.trapz(self.phi, self.g.x1_grid)
# 
#     def work(self):
#         u = self.g.x3(self.g.x1_grid)
#         if self.load_type == 'concentrated':
#             self.W = self.load*u[-1]
#         elif self.load_type == 'distributed':
#             self.W = np.trapz(self.load*u, self.g.x1_grid)
# 
#     def residual(self):
#         self.R = self.U - self.W
# 
#     def minimum_potential(self, x0=[0,0,0]):
#         def to_optimize(x):
#             self.g.D = [0,0] + list(x)
#             self.work()
#             self.free_energy()
#             self.strain_energy()
#             self.residual()
#             return self.R
# 
#         # With bounds
#         bounds = np.array(((-0.01,0.01),)*len(x0))
# 
#         res = minimize(to_optimize, x0, bounds=bounds)
#         self.g.D = [0,0] + list(res.x)
#         self.work()
#         self.free_energy()
#         self.strain_energy()
#         self.residual()
#         return(res.x, res.fun)


# Path: examples/structural/beam_theory/beam_debugging.py
from aeropy.geometry.parametric import CoordinateSystem
from aeropy.structural.stable_solution import (properties, boundary_conditions, euler_bernoulle)
from aeropy.structural.shell import shell
from aeropy.xfoil_module import output_reader
from scipy.integrate import quad
from scipy.optimize import minimize
from optimization_tools.DOE import DOE
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np
import pickle



# beam.minimum_potential(x0=[0,0], input_function = lambda x: [0,0] + list(x),
                       # bounds = bounds)

# Euler Bernoulle
EB = euler_bernoulle(bp, -1, 'concentrated', beam.g_c, x=beam.g_c.x1_grid)
EB.analytical_solutions()
EB.bending_strain()
EB.free_energy()
EB.strain_energy()
print('U_EB: ', EB.U)
print('U_PL: ', beam.U)
# mod = 1/(1-0.3**2)
# print('U_corrected: ', beam.U/mod)
#
print('Bending_strain_EB: ', EB.B)
print('Bending_strain_PL: ', beam.rho)
#
# print('Phi_EB: ', EB.phi)
# print('Phi_PL: ', beam.phi)
# print('Phi_PL_B: ', beam.phi_B)
#
# print('ddx3_EB: ', EB.g.x3(beam.g_c.x1_grid, 'x11'))
# print('ddx3_PL: ', beam.g_c.x3(beam.g_c.x1_grid, 'x11'))

print('theta1',beam.g_c.x1(beam.g_c.x1_grid, 'theta1'))
print('theta11', beam.g_c.x1(beam.g_c.x1_grid, 'theta11'))
P = -1
E = bp.young

L = 1