<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 For now, neglect rotational inertia.
Interpolation functions
Step1:  Mass matrix
Define the density distribution (linear)
Step2:  Integrate the density distribution with the shape functions.
Step3:  Special case
Step4:  Shape integrals
As well as the actual mass matrix, the shape integrals are needed for the multibody dynamics equations
Step5:  First shape integral
Step6:  The mass matrix should be the same as the trace of the 2nd shape integrals
Step7:  Stiffness matrix
Differentiate the shape functions
Step8:  Define the stiffness distribution (linear)
Step9:  Note that $EI_y$ refers to the stiffness for bending in the $y$-direction, not about the $y$ axis.
Try and simplify results by using the average values where they come up
Step10:  Special case
Step11:  This is the same as Rao2004 p.326
Stress stiffening
Now derive the stress stiffening (e.g. centrifugal stiffening) matrix. Ref Cook1989.
Need the matrix $\boldsymbol{G}$ such that
$$
\begin{bmatrix} u_{,x} \ v_{,x} \ w_{,x} \ \end{bmatrix} = \boldsymbol{Gq}
$$
Step12:  This agrees with Cook1989, p.434, for the transverse directions. 
Generalised forces
Define the force distribution (linear)
Step13:  First shape integral for applied forces
Step14:  Second shape integral for applied forces
Step15:  The generalised nodal forces are given by the trace of this (the other parts of it are used to find the moments on the whole body directly...)
Step16:  Special case -- constant force
Step17:  Another example -- a uniform distributed force in the z direction
Step19:  Output
Write a Python file with functions to calculate all of these
<ASSISTANT_TASK:>
Python Code:
xi, l, rho = symbols('xi, l, rho')
# Shape functions
S = Matrix(np.zeros((4, 12)))
x2 = (1 - xi)
S[0, 0 ] =  x2                     # extension
S[0, 6 ] =  xi
S[1, 1 ] =  x2**2 * (3 - 2*x2)     # y-deflection
S[1, 7 ] =  xi**2 * (3 - 2*xi)
S[1, 5 ] = -x2**2 * (x2 - 1) * l  
S[1, 11] =  xi**2 * (xi - 1) * l  
S[2, 2 ] =  x2**2 * (3 - 2*x2)     # z-deflection
S[2, 8 ] =  xi**2 * (3 - 2*xi)
S[2, 4 ] =  x2**2 * (x2 - 1) * l
S[2, 10] = -xi**2 * (xi - 1) * l
S[3, 3 ] =  x2                     # torsion
S[3, 9 ] =  xi
#S[4, 2 ] =  6 * x2 * (x2 - 1) / l  # y-rotation
#S[4, 8 ] =  6 * xi * (xi - 1) / l
#S[4, 4 ] = -x2 * (3*x2 - 2)  
#S[4, 10] =  xi * (3*xi - 2)  
#S[5, 1 ] = -6 * x2 * (x2 - 1) / l  # z-rotation
#S[5, 7 ] = -6 * xi * (xi - 1) / l
#S[5, 5 ] =  x2 * (3*x2 - 2)  
#S[5, 11] =  xi * (3*xi - 2)  
S[:3, :].T
titles = ['x-defl', 'y-defl', 'z-defl', 'torsion']
for i in range(4):
    sympy.plot(*([xx.subs(l, 2) for xx in S[i,:] if xx != 0] + [(xi, 0, 1)]), 
               title=titles[i])
Explanation: For now, neglect rotational inertia.
Interpolation functions
End of explanation
rho1, rho2 = symbols('rho_1, rho_2')
rho = (1 - xi)*rho1 + xi*rho2
rho
Explanation: Mass matrix
Define the density distribution (linear):
End of explanation
def sym_me():
    m = Matrix(np.diag([rho, rho, rho, 0]))
    integrand = S.T * m * S
    me = integrand.applyfunc(
        lambda xxx: l * sympy.integrate(xxx, (xi, 0, 1)).expand().factor()
    )
    return me
me = sym_me()
me.shape
me[0,:]
me[6,:]
me[1,:]
Explanation: Integrate the density distribution with the shape functions.
End of explanation
me.subs({rho2: rho1})/rho1
Explanation: Special case: rho1 == rho2
End of explanation
mass = l * sympy.integrate(rho, (xi, 0, 1)).factor()
mass
Explanation: Shape integrals
As well as the actual mass matrix, the shape integrals are needed for the multibody dynamics equations:
\begin{align}
m                   &= \int \mathrm{d}m      \
\boldsymbol{S}      &= \int \boldsymbol{S} \mathrm{d}m  \
\boldsymbol{S}_{kl} &= \int \boldsymbol{S}_k^T \boldsymbol{S}_l \mathrm{d}m
\end{align}
where $\boldsymbol{S}_k$ is the $k$th row of the element shape function.
The mass is the average density times the length:
End of explanation
shape_integral_1 = S[:3, :].applyfunc(
    lambda xxx: l * sympy.integrate(rho * xxx, (xi, 0, 1)).expand().simplify()
)
shape_integral_1.T
shape_integral_2 = [
    [l * (S[i, :].T * S[j, :]).applyfunc(
        lambda xxx: sympy.integrate(rho * xxx, (xi, 0, 1)).expand().simplify())
     for j in range(3)]
    for i in range(3)
]
Explanation: First shape integral:
End of explanation
(shape_integral_2[0][0] + shape_integral_2[1][1] + shape_integral_2[2][2] - me).expand()
Explanation: The mass matrix should be the same as the trace of the 2nd shape integrals:
End of explanation
B = Matrix(np.zeros((4, 12)))
B[0, :] = S[0, :].diff(xi, 1) / l
B[1, :] = S[1, :].diff(xi, 2) / l**2
B[2, :] = S[2, :].diff(xi, 2) / l**2
B[3, :] = S[3, :].diff(xi, 1) / l
B.simplify()
B[:, :6]
B[:, 6:]
titles = ['x-defl', 'y-defl', 'z-defl', 'torsion']
for i in range(4):
    sympy.plot(*([xx.subs(l, 2) for xx in B[i,:] if xx != 0] + [(xi, 0, 1)]), 
               title=titles[i])
Explanation: Stiffness matrix
Differentiate the shape functions:
End of explanation
EA1, EA2, EIy1, EIy2, EIz1, EIz2, GJ1, GJ2 = symbols('EA_1, EA_2, EIy_1, EIy_2, EIz_1, EIz_2, GJ_1, GJ_2')
EA  = (1 - xi)*EA1  + xi*EA2
EIy = (1 - xi)*EIy1 + xi*EIy2
EIz = (1 - xi)*EIz1 + xi*EIz2
GJ  = (1 - xi)*GJ1  + xi*GJ2
EA, EIy, EIz, GJ
Explanation: Define the stiffness distribution (linear):
End of explanation
Ex, Ey, Ez, Gx = symbols('E_x, E_y, E_z, G_x')
def sym_ke():
    # Note the order -- y deflections depend on EIz etc
    E = Matrix(np.diag([EA, EIz, EIy, GJ]))
    integrand = B.T * E * B
    ke = integrand.applyfunc(
        lambda xxx: l * sympy.integrate(xxx, (xi, 0, 1)).factor() #.subs((EA1+EA2), 2*Ex) #.expand().factor()
    )
    return ke
ke = sym_ke()
def simplify_ke(ke, EI=False):
    result = ke.applyfunc(
        lambda xxx: xxx.subs((EA1+EA2), 2*Ex).subs((GJ1+GJ2), 2*Gx))
    if EI:
        result = result.applyfunc(
            lambda xxx: xxx.subs((EIy1+EIy2), 2*Ey).subs((EIz1+EIz2), 2*Ez))
    return result
kem = simplify_ke(ke, True)
kem[:, :6]
kem[:, 6:]
Explanation: Note that $EI_y$ refers to the stiffness for bending in the $y$-direction, not about the $y$ axis.
Try and simplify results by using the average values where they come up
End of explanation
kem.subs({EA1: Ex, EA2: Ex, EIy1: Ey, EIy2: Ey, EIz1: Ez, EIz2: Ez})
Explanation: Special case: uniform stiffness.
End of explanation
# G = Matrix(np.zeros((3, 12)))
# G[0, :] = S[0, :].diff(xi, 1) / l
# G[1, :] = S[1, :].diff(xi, 1) / l
# G[2, :] = S[2, :].diff(xi, 1) / l
G = S[:3, :].diff(xi, 1) / l
G.simplify()
G[:, :6]
G[:, 6:]
def sym_ks_axial_force():
    # Unit axial force (absorbing area from integral), block matrix 3 times
    #smat3 = Matrix(9, 9, lambda i, j: 1 if i == j and i % 3 == 0 else 0)
    #integrand = G.T * smat3 * G
    integrand = G.T * G
    ks = integrand.applyfunc(
        lambda xxx: l * sympy.integrate(xxx, (xi, 0, 1)).factor() #.subs((EA1+EA2), 2*Ex) #.expand().factor()
    )
    return ks
ks = sym_ks_axial_force()
ks * 30*l
Explanation: This is the same as Rao2004 p.326
Stress stiffening
Now derive the stress stiffening (e.g. centrifugal stiffening) matrix. Ref Cook1989.
Need the matrix $\boldsymbol{G}$ such that
$$
\begin{bmatrix} u_{,x} \ v_{,x} \ w_{,x} \ \end{bmatrix} = \boldsymbol{Gq}
$$
End of explanation
fx1, fx2, fy1, fy2, fz1, fz2 = symbols('f_x1, f_x2, f_y1, f_y2, f_z1, f_z2')
fx = (1 - xi)*fx1 + xi*fx2
fy = (1 - xi)*fy1 + xi*fy2
fz = (1 - xi)*fz1 + xi*fz2
f = Matrix([fx, fy, fz])
f
# Shape functions for applied force -- linear
SF = Matrix(np.zeros((3, 12)))
SF[0, 0 ] =  x2                     # x
SF[0, 6 ] =  xi
SF[1, 1 ] =  x2                     # y
SF[1, 7 ] =  xi
SF[2, 2 ] =  x2                     # z
SF[2, 8 ] =  xi
SF
Explanation: This agrees with Cook1989, p.434, for the transverse directions. 
Generalised forces
Define the force distribution (linear):
End of explanation
shape_integral_F1 = SF.applyfunc(
    lambda xxx: l * sympy.integrate(xxx, (xi, 0, 1)).expand().simplify()
)
shape_integral_F1
Explanation: First shape integral for applied forces:
End of explanation
shape_integral_F2 = [
    [l * (S[i, :].T * SF[j, :]).applyfunc(
        lambda xxx: sympy.integrate(xxx, (xi, 0, 1)).expand().simplify())
     for j in range(3)]
    for i in range(3)
]
Explanation: Second shape integral for applied forces:
End of explanation
F = shape_integral_F2[0][0] + shape_integral_F2[1][1] + shape_integral_F2[2][2]
Explanation: The generalised nodal forces are given by the trace of this (the other parts of it are used to find the moments on the whole body directly...):
End of explanation
F * Matrix([fx1, fy1, fz1, 0, 0, 0, fx1, fy1, fz1, 0, 0, 0])
Explanation: Special case -- constant force:
End of explanation
F * Matrix([0, 0, fz1, 0, 0, 0, 0, 0, fz1, 0, 0, 0]) / (l/12*fz1)
Explanation: Another example -- a uniform distributed force in the z direction:
End of explanation
def numpy_array_str(expr):
    return str(expr) \
        .replace('Matrix([', 'array([\n') \
        .replace('], [', '],\n[') \
        .replace(']])', ']\n])')
        
def numpy_array_str_2x2(arr):
    return ',\n'.join(['[{}]'.format(',\n'.join([numpy_array_str(arr[i][j])
                                                 for j in range(3)]))
                        for i in range(3)])
import datetime
code = 
# Automatically generated from SymPy in notebook
# {date}
from __future__ import division
from numpy import array
def mass(l, rho_1, rho_2):
    return {mass}
    
def S1(l, rho_1, rho_2):
    return {S1}
    
def S2(l, rho_1, rho_2):
    return [
        {S2}
    ]
    
def K(l, E_x, G_x, EIy_1, EIy_2, EIz_1, EIz_2):
    return {K}
    
def Ks(l):
    return {Ks}
    
def F1(l):
    return {F1}
    
def F2(l):
    return [
        {F2}
    ]
    
.format(
    date=datetime.datetime.now(),
    mass=mass,
    S1=numpy_array_str(shape_integral_1),
    S2=numpy_array_str_2x2(shape_integral_2),
    K=numpy_array_str(simplify_ke(ke)),
    Ks=numpy_array_str(ks),
    F1=numpy_array_str(shape_integral_F1),
    F2=numpy_array_str_2x2(shape_integral_F2)
)
with open('../beamfe/tapered_beam_element_integrals.py', 'wt') as f:
    f.write(code)
Explanation: Output
Write a Python file with functions to calculate all of these:
End of explanation
<END_TASK>