# Source code for openmdao.test_suite.components.three_bar_truss

"""
3 Bar Truss Problem.

- 3 continuous design variables -  Cross sectional area of each bar.

- 3 discrete material choices for each bar. Complete enumeration requires 64 continuous optimization.

The global optima is 3,3,x (x can be anything), with mass = 5.28 kg
"""

import numpy as np

from openmdao.core.explicitcomponent import ExplicitComponent

[docs]def stress_calc(A, E):
"""
Calculate stress in a 3-bar truss given area and Young's modulus.

Parameters
----------
A : float or ndarray
Area of the bar.
E : float or ndarray
Young's Modules of the bar.

Returns
-------
float or ndarray
Stress in the bar.
"""
P  = 120000.0
L = np.array([np.sqrt(1.2**2 + 1.2**2), 1.2, np.sqrt(1.2**2 + 1.2**2)])

#Local stiffness matrix
theta1 = -45.0*np.pi/180.0
theta2 = -90.0*np.pi/180.0
theta3 = -135*np.pi/180.0

K0 = (E[0]*A[0]/L[0])*np.dot(np.array([[np.cos(theta1), np.sin(theta1)]]).T,
np.array([[np.cos(theta1), np.sin(theta1)]]))
K1 = (E[1]*A[1]/L[1])*np.dot(np.array([[np.cos(theta2), np.sin(theta2)]]).T,
np.array([[np.cos(theta2), np.sin(theta2)]]))
K2 = (E[2]*A[2]/L[2])*np.dot(np.array([[np.cos(theta3), np.sin(theta3)]]).T,
np.array([[np.cos(theta3), np.sin(theta3)]]))

# Global (total) stiffness matrix
K = K0 + K1 + K2

theta4 = -65.0*np.pi/180.0
p = P*np.array([[np.cos(theta4), np.sin(theta4)]]).T
# Displacement matrix
u = np.dot(np.linalg.inv(K), p)

#Delta change in length
DL = np.zeros([3])
DL[0] = np.sqrt((-L[0]*np.cos(theta1) - u[0])**2 + (-L[0]*np.sin(theta1) - u[1])**2) - L[0]
DL[1] = np.sqrt((-L[1]*np.cos(theta2) - u[0])**2 + (-L[1]*np.sin(theta2) - u[1])**2) - L[1]
DL[2] = np.sqrt((-L[2]*np.cos(theta3) - u[0])**2 + (-L[2]*np.sin(theta3) - u[1])**2) - L[2]

#Stress in each element
sigma = E*DL/L
return sigma

[docs]class ThreeBarTruss(ExplicitComponent):
"""
3 Bar truss problem with 3 continuous design variables and 3 discrete
material choices. Material chosen as follows:

1 Aluminum
2 Titanium
3 Steel
4 Nickel

This component calculates the total mass of the truss.
"""

[docs]    def setup(self):
# Continuous Inputs
desc='Cross-sectional area of beam 1')
desc='Cross-sectional area of beam 2')
desc='Cross-sectional area of beam 3')

# Discrete Inputs
desc='Material ID of beam 1')
desc='Material ID of beam 2')
desc='Material ID of beam 3')

# Outputs

self.rho = { 1 : 2700.0,
2 : 4500.0,
3 : 7872.0,
4 : 8800.0 }

self.E = { 1 : 68.9e9,
2 : 116.0e9,
3 : 205.0e9,
4 : 207.0e9 }

self.sigma_y = { 1 : 55.2e6,
2 : 140.0e6,
3 : 285.0e6,
4 : 59.0e6 }

self.declare_partials(of='*', wrt='area*', method='fd')

[docs]    def compute(self, inputs, outputs):
"""
Define the function f(xI, xC).

Here xI is integer and xC is continuous.
"""

# Convert areas to m**2
area1 = inputs['area1']*.0001
area2 = inputs['area2']*.0001
area3 = inputs['area3']*.0001
mat1 = int(inputs['mat1'])
mat2 = int(inputs['mat2'])
mat3 = int(inputs['mat3'])
area1 *= 2.

len1 = np.sqrt(1.2**2 + 1.2**2)
len2 = 1.2
len3 = np.sqrt(1.2**2 + 1.2**2)

rho1 = self.rho[mat1]
rho2 = self.rho[mat2]
rho3 = self.rho[mat3]

outputs['mass'] = rho1*area1*len1 + rho2*area2*len2 + rho3*area3*len3

E1 = self.E[mat1]
E2 = self.E[mat2]
E3 = self.E[mat3]
sigma_y1 = self.sigma_y[mat1]
sigma_y2 = self.sigma_y[mat2]
sigma_y3 = self.sigma_y[mat3]

E = np.array([E1, E2, E3])
A = np.array([area1, area2, area3])
sigma_y = np.array([sigma_y1, sigma_y2, sigma_y3])

sigma = stress_calc(A, E)
outputs['stress'] = (np.abs(sigma)/sigma_y)

[docs]class ThreeBarTrussVector(ExplicitComponent):
"""
3 Bar truss problem with 3 continuous design variables and 3 discrete
material choices. Material chosen as follows:

1 Aluminum
2 Titanium
3 Steel
4 Nickel

This component calculates the total mass of the truss.

This version places all areas in a single vector and all materials in
a single vector for test purporses.
"""

[docs]    def setup(self):
# Continuous Inputs
desc='Cross-sectional area of beams')

# Discrete Inputs
self.add_input('mat', np.array([1, 1, 1]), #lower=1, upper=4,
desc='Material ID of beams')

# Outputs

self.rho = { 1 : 2700.0,
2 : 4500.0,
3 : 7872.0,
4 : 8800.0 }

self.E = { 1 : 68.9e9,
2 : 116.0e9,
3 : 205.0e9,
4 : 207.0e9 }

self.sigma_y = { 1 : 55.2e6,
2 : 140.0e6,
3 : 285.0e6,
4 : 59.0e6 }

self.declare_partials(of='*', wrt='area', method='fd')

[docs]    def compute(self, inputs, outputs):
""" Define the function f(xI, xC)
Here xI is integer and xC is continuous"""

# Area convert to m**2
area = inputs['area']*.0001
mat = inputs['mat']

length = np.array([np.sqrt(1.2**2 + 1.2**2),
1.2,
np.sqrt(1.2**2 + 1.2**2)])

rho = np.zeros((3, ))
try:
rho[0] = self.rho[mat[0]]
except:
pass
rho[1] = self.rho[mat[1]]
rho[2] = self.rho[mat[2]]

outputs['mass'] = np.sum(rho*area*length)

E = np.zeros((3, ))
E[0] = self.E[mat[0]]
E[1] = self.E[mat[1]]
E[2] = self.E[mat[2]]
sigma_y = np.zeros((3, ))

sigma_y[0] = self.sigma_y[mat[0]]
sigma_y[1] = self.sigma_y[mat[1]]
sigma_y[2] = self.sigma_y[mat[2]]

sigma = stress_calc(area, E)
outputs['stress'] = (np.abs(sigma)/sigma_y)