from sympy import *
from sympy.vector import CoordSys3D
N = CoordSys3D('N')
x1, x2, x3 = symbols("x_1 x_2 x_3")
alpha1, alpha2, alpha3 = symbols("alpha_1 alpha_2 alpha_3")
R, L, ga, gv = symbols("R L g_a g_v")
init_printing()

"""
Explanation: Corrugated Shells
Init symbols for sympy
End of explanation
"""


a1 = pi / 2 + (L / 2 - alpha1)/R

x = (R + alpha3 + ga * cos(gv * a1)) * cos(a1)
y = alpha2
z = (R + alpha3 + ga * cos(gv * a1)) * sin(a1)

r = x*N.i + y*N.j + z*N.k

"""
Explanation: Corrugated cylindrical coordinates
End of explanation
"""


R1=r.diff(alpha1)
R2=r.diff(alpha2)
R3=r.diff(alpha3)

trigsimp(R1)

R2

R3

"""
Explanation: Base Vectors $\vec{R}_1, \vec{R}_2, \vec{R}_3$
End of explanation
"""


eps=trigsimp(R1.dot(R2.cross(R3)))
R_1=simplify(trigsimp(R2.cross(R3)/eps))
R_2=simplify(trigsimp(R3.cross(R1)/eps))
R_3=simplify(trigsimp(R1.cross(R2)/eps))

R_1

R_2

R_3

"""
Explanation: Base Vectors $\vec{R}^1, \vec{R}^2, \vec{R}^3$
End of explanation
"""


dx1da1=R1.dot(N.i)
dx1da2=R2.dot(N.i)
dx1da3=R3.dot(N.i)

dx2da1=R1.dot(N.j)
dx2da2=R2.dot(N.j)
dx2da3=R3.dot(N.j)

dx3da1=R1.dot(N.k)
dx3da2=R2.dot(N.k)
dx3da3=R3.dot(N.k)

A=Matrix([[dx1da1, dx1da2, dx1da3], [dx2da1, dx2da2, dx2da3], [dx3da1, dx3da2, dx3da3]])
simplify(A)

w=Function('w')
z=Function('z')
a=Function('a')

dx1da1 = -z(alpha1)*sin((L/2-alpha1)/R)+w(alpha1, alpha3)*cos((L/2-alpha1)/R)
dx1da2 = 0
dx1da3 = -sin((L/2-alpha1)/R)

dx2da1 = 0
dx2da2 = 1
dx2da3 = 0

dx3da1 = z(alpha1)*cos((L/2-alpha1)/R)+w(alpha1, alpha3)*sin((L/2-alpha1)/R)
dx3da2 = 0
dx3da3 = cos((L/2-alpha1)/R)

A_s=Matrix([[dx1da1, dx1da2, dx1da3], [dx2da1, dx2da2, dx2da3], [dx3da1, dx3da2, dx3da3]])
A_s

A_inv = A_s**-1
A_inv=simplify(A_inv)
A_inv

trigsimp(A_s.det())

"""
Explanation: Jacobi matrix:
$ A = \left( 
\begin{array}{ccc} 
\frac{\partial x_1}{\partial \alpha_1} & \frac{\partial x_1}{\partial \alpha_2} & \frac{\partial x_1}{\partial \alpha_3} \
\frac{\partial x_2}{\partial \alpha_1} & \frac{\partial x_2}{\partial \alpha_2} & \frac{\partial x_3}{\partial \alpha_3} \
\frac{\partial x_3}{\partial \alpha_1} & \frac{\partial x_3}{\partial \alpha_2} & \frac{\partial x_3}{\partial \alpha_3} \
\end{array} 
\right)$
$ \left[
\begin{array}{ccc} 
\vec{R}_1 & \vec{R}_2 & \vec{R}_3
\end{array} 
\right] = \left[
\begin{array}{ccc} 
\vec{e}_1 & \vec{e}_2 & \vec{e}_3
\end{array} 
\right] \cdot  \left( 
\begin{array}{ccc} 
\frac{\partial x_1}{\partial \alpha_1} & \frac{\partial x_1}{\partial \alpha_2} & \frac{\partial x_1}{\partial \alpha_3} \
\frac{\partial x_2}{\partial \alpha_1} & \frac{\partial x_2}{\partial \alpha_2} & \frac{\partial x_3}{\partial \alpha_3} \
\frac{\partial x_3}{\partial \alpha_1} & \frac{\partial x_3}{\partial \alpha_2} & \frac{\partial x_3}{\partial \alpha_3} \
\end{array} 
\right) = \left[
\begin{array}{ccc} 
\vec{e}_1 & \vec{e}_2 & \vec{e}_3
\end{array} 
\right] \cdot A$
$ \left[
\begin{array}{ccc} 
\vec{e}_1 & \vec{e}_2 & \vec{e}_3
\end{array} 
\right] =\left[
\begin{array}{ccc} 
\vec{R}_1 & \vec{R}_2 & \vec{R}_3
\end{array} 
\right] \cdot A^{-1}$
End of explanation
"""


g11=R1.dot(R1)
g12=R1.dot(R2)
g13=R1.dot(R3)

g21=R2.dot(R1)
g22=R2.dot(R2)
g23=R2.dot(R3)

g31=R3.dot(R1)
g32=R3.dot(R2)
g33=R3.dot(R3)

G=Matrix([[g11, g12, g13],[g21, g22, g23], [g31, g32, g33]])
G=trigsimp(G)
G

G_s=A_s.T*A_s
simplify(G_s)

"""
Explanation: Metric tensor
${\displaystyle \hat{G}=\sum_{i,j} g^{ij}\vec{R}_i\vec{R}_j}$
End of explanation
"""


g_11=R_1.dot(R_1)
g_12=R_1.dot(R_2)
g_13=R_1.dot(R_3)

g_21=R_2.dot(R_1)
g_22=R_2.dot(R_2)
g_23=R_2.dot(R_3)

g_31=R_3.dot(R_1)
g_32=R_3.dot(R_2)
g_33=R_3.dot(R_3)

G_con=Matrix([[g_11, g_12, g_13],[g_21, g_22, g_23], [g_31, g_32, g_33]])
G_con=trigsimp(G_con)
G_con

G_con_s=G_s**-1
simplify(G_con_s)

"""
Explanation: ${\displaystyle \hat{G}=\sum_{i,j} g_{ij}\vec{R}^i\vec{R}^j}$
End of explanation
"""


dR1dalpha1 = trigsimp(R1.diff(alpha1))
dR1dalpha1

"""
Explanation: Derivatives of vectors
Derivative of base vectors
End of explanation
"""


dR1dalpha2 = trigsimp(R1.diff(alpha2))
dR1dalpha2

dR1dalpha3 = trigsimp(R1.diff(alpha3))
dR1dalpha3

"""
Explanation: $ \frac { d\vec{R_1} } { d\alpha_1} = -\frac {1}{R} \left( 1+\frac{\alpha_3}{R} \right) \vec{R_3} $
End of explanation
"""


dR2dalpha1 = trigsimp(R2.diff(alpha1))
dR2dalpha1

dR2dalpha2 = trigsimp(R2.diff(alpha2))
dR2dalpha2

dR2dalpha3 = trigsimp(R2.diff(alpha3))
dR2dalpha3

dR3dalpha1 = trigsimp(R3.diff(alpha1))
dR3dalpha1

"""
Explanation: $ \frac { d\vec{R_1} } { d\alpha_3} = \frac {1}{R} \frac {1}{1+\frac{\alpha_3}{R}}  \vec{R_1} $
End of explanation
"""


dR3dalpha2 = trigsimp(R3.diff(alpha2))
dR3dalpha2

dR3dalpha3 = trigsimp(R3.diff(alpha3))
dR3dalpha3

"""
Explanation: $ \frac { d\vec{R_3} } { d\alpha_1} = \frac {1}{R} \frac {1}{1+\frac{\alpha_3}{R}}  \vec{R_1} $
End of explanation
"""


u1=Function('u^1')
u2=Function('u^2')
u3=Function('u^3')
q=Function('q') # q(alpha3) = 1+alpha3/R
K = Symbol('K') # K = 1/R

u1_nabla1 = u1(alpha1, alpha2, alpha3).diff(alpha1) + u3(alpha1, alpha2, alpha3) * K / q(alpha3)
u2_nabla1 = u2(alpha1, alpha2, alpha3).diff(alpha1)
u3_nabla1 = u3(alpha1, alpha2, alpha3).diff(alpha1) - u1(alpha1, alpha2, alpha3) * K * q(alpha3)

u1_nabla2 = u1(alpha1, alpha2, alpha3).diff(alpha2)
u2_nabla2 = u2(alpha1, alpha2, alpha3).diff(alpha2)
u3_nabla2 = u3(alpha1, alpha2, alpha3).diff(alpha2)

u1_nabla3 = u1(alpha1, alpha2, alpha3).diff(alpha3) + u1(alpha1, alpha2, alpha3) * K / q(alpha3)
u2_nabla3 = u2(alpha1, alpha2, alpha3).diff(alpha3)
u3_nabla3 = u3(alpha1, alpha2, alpha3).diff(alpha3)
# $\nabla_2 u^2 = \frac { \partial u^2 } { \partial \alpha_2}$

grad_u = Matrix([[u1_nabla1, u2_nabla1, u3_nabla1],[u1_nabla2, u2_nabla2, u3_nabla2], [u1_nabla3, u2_nabla3, u3_nabla3]])
grad_u

G_s = Matrix([[q(alpha3)**2, 0, 0],[0, 1, 0], [0, 0, 1]])
grad_u_down=grad_u*G_s
expand(simplify(grad_u_down))

"""
Explanation: $ \frac { d\vec{R_3} } { d\alpha_3} = \vec{0} $
Derivative of vectors
$ \vec{u} = u^1 \vec{R_1} + u^2\vec{R_2} + u^3\vec{R_3} $
$ \frac { d\vec{u} } { d\alpha_1} =  \frac { d(u^1\vec{R_1}) } { d\alpha_1} + \frac { d(u^2\vec{R_2}) } { d\alpha_1}+ \frac { d(u^3\vec{R_3}) } { d\alpha_1} = \frac { du^1 } { d\alpha_1} \vec{R_1} + u^1 \frac { d\vec{R_1} } { d\alpha_1} + \frac { du^2 } { d\alpha_1} \vec{R_2} + u^2 \frac { d\vec{R_2} } { d\alpha_1} + \frac { du^3 } { d\alpha_1} \vec{R_3} + u^3 \frac { d\vec{R_3} } { d\alpha_1} = \frac { du^1 } { d\alpha_1} \vec{R_1} - u^1 \frac {1}{R} \left( 1+\frac{\alpha_3}{R} \right) \vec{R_3} + \frac { du^2 } { d\alpha_1} \vec{R_2}+ \frac { du^3 } { d\alpha_1} \vec{R_3} + u^3 \frac {1}{R} \frac {1}{1+\frac{\alpha_3}{R}}  \vec{R_1}$
Then
$ \frac { d\vec{u} } { d\alpha_1} = \left( \frac { du^1 } { d\alpha_1} + u^3 \frac {1}{R} \frac {1}{1+\frac{\alpha_3}{R}}  \right) \vec{R_1} + \frac { du^2 } { d\alpha_1} \vec{R_2} + \left( \frac { du^3 } { d\alpha_1} - u^1 \frac {1}{R} \left( 1+\frac{\alpha_3}{R} \right) \right) \vec{R_3}$
$ \frac { d\vec{u} } { d\alpha_2} =  \frac { d(u^1\vec{R_1}) } { d\alpha_2} + \frac { d(u^2\vec{R_2}) } { d\alpha_2}+ \frac { d(u^3\vec{R_3}) } { d\alpha_2} =  \frac { du^1 } { d\alpha_2} \vec{R_1} + \frac { du^2 } { d\alpha_2} \vec{R_2} + \frac { du^3 } { d\alpha_2} \vec{R_3} $
Then
$ \frac { d\vec{u} } { d\alpha_2} = \frac { du^1 } { d\alpha_2} \vec{R_1} + \frac { du^2 } { d\alpha_2} \vec{R_2} + \frac { du^3 } { d\alpha_2} \vec{R_3} $
$ \frac { d\vec{u} } { d\alpha_3} = \frac { d(u^1\vec{R_1}) } { d\alpha_3} + \frac { d(u^2\vec{R_2}) } { d\alpha_3}+ \frac { d(u^3\vec{R_3}) } { d\alpha_3} = 
 \frac { du^1 } { d\alpha_3} \vec{R_1} + u^1 \frac { d\vec{R_1} } { d\alpha_3} + \frac { du^2 } { d\alpha_3} \vec{R_2} + u^2 \frac { d\vec{R_2} } { d\alpha_3} + \frac { du^3 } { d\alpha_3} \vec{R_3} + u^3 \frac { d\vec{R_3} } { d\alpha_3} = \frac { du^1 } { d\alpha_3} \vec{R_1} + u^1 \frac {1}{R} \frac {1}{1+\frac{\alpha_3}{R}}  \vec{R_1} + \frac { du^2 } { d\alpha_3} \vec{R_2}+ \frac { du^3 } { d\alpha_3} \vec{R_3} $
Then
$ \frac { d\vec{u} } { d\alpha_3} = \left( \frac { du^1 } { d\alpha_3}  + u^1 \frac {1}{R} \frac {1}{1+\frac{\alpha_3}{R}} \right) \vec{R_1} + \frac { du^2 } { d\alpha_3} \vec{R_2}+ \frac { du^3 } { d\alpha_3} \vec{R_3}$
Gradient of vector
$\nabla_1 u^1 = \frac { \partial u^1 } { \partial \alpha_1} + u^3 \frac {1}{R} \frac {1}{1+\frac{\alpha_3}{R}}$
$\nabla_1 u^2 = \frac { \partial u^2 } { \partial \alpha_1} $
$\nabla_1 u^3 = \frac { \partial u^3 } { \partial \alpha_1} - u^1 \frac {1}{R} \left( 1+\frac{\alpha_3}{R} \right) $
$\nabla_2 u^1 = \frac { \partial u^1 } { \partial \alpha_2}$
$\nabla_2 u^2 = \frac { \partial u^2 } { \partial \alpha_2}$
$\nabla_2 u^3 = \frac { \partial u^3 } { \partial \alpha_2}$
$\nabla_3 u^1 = \frac { \partial u^1 } { \partial \alpha_3} + u^1 \frac {1}{R} \frac {1}{1+\frac{\alpha_3}{R}}$
$\nabla_3 u^2 = \frac { \partial u^2 } { \partial \alpha_3} $
$\nabla_3 u^3 = \frac { \partial u^3 } { \partial \alpha_3}$
$ \nabla \vec{u} = \left( 
\begin{array}{ccc} 
\nabla_1 u^1 & \nabla_1 u^2 & \nabla_1 u^3 \
\nabla_2 u^1 & \nabla_2 u^2 & \nabla_2 u^3 \
\nabla_3 u^1 & \nabla_3 u^2 & \nabla_3 u^3 \
\end{array} 
\right)$
End of explanation
"""


B = zeros(9, 12)
B[0,1] = (1+alpha3/R)**2
B[0,8] = (1+alpha3/R)/R

B[1,2] = (1+alpha3/R)**2

B[2,0] = (1+alpha3/R)/R
B[2,3] = (1+alpha3/R)**2

B[3,5] = S(1)
B[4,6] = S(1)
B[5,7] = S(1)

B[6,9] = S(1)
B[6,0] = -(1+alpha3/R)/R
B[7,10] = S(1)
B[8,11] = S(1)

B

"""
Explanation: $ 
\left( 
\begin{array}{c} 
\nabla_1 u_1 \ \nabla_2 u_1 \ \nabla_3 u_1 \
\nabla_1 u_2 \ \nabla_2 u_2 \ \nabla_3 u_2 \
\nabla_1 u_3 \ \nabla_2 u_3 \ \nabla_3 u_3 \
\end{array} 
\right) 
=
\left( 
\begin{array}{c}
\left( 1+\frac{\alpha_2}{R} \right)^2  \frac { \partial u^1 } { \partial \alpha_1} + u^3 \frac {\left( 1+\frac{\alpha_3}{R} \right)}{R} \
\left( 1+\frac{\alpha_2}{R} \right)^2 \frac { \partial u^1 } { \partial \alpha_2} \
\left( 1+\frac{\alpha_3}{R} \right)^2  \frac { \partial u^1 } { \partial \alpha_3} + u^1 \frac {\left( 1+\frac{\alpha_3}{R} \right)}{R} \
\frac { \partial u^2 } { \partial \alpha_1} \
\frac { \partial u^2 } { \partial \alpha_2} \
\frac { \partial u^2 } { \partial \alpha_3} \
\frac { \partial u^3 } { \partial \alpha_1} - u^1 \frac {\left( 1+\frac{\alpha_3}{R} \right)}{R} \
\frac { \partial u^3 } { \partial \alpha_2} \
\frac { \partial u^3 } { \partial \alpha_3} \
\end{array} 
\right)
$
$ 
\left( 
\begin{array}{c} 
\nabla_1 u_1 \ \nabla_2 u_1 \ \nabla_3 u_1 \
\nabla_1 u_2 \ \nabla_2 u_2 \ \nabla_3 u_2 \
\nabla_1 u_3 \ \nabla_2 u_3 \ \nabla_3 u_3 \
\end{array} 
\right)
=
B \cdot
\left( 
\begin{array}{c} 
u^1 \
\frac { \partial u^1 } { \partial \alpha_1} \
\frac { \partial u^1 } { \partial \alpha_2} \
\frac { \partial u^1 } { \partial \alpha_3} \
u^2 \
\frac { \partial u^2 } { \partial \alpha_1} \
\frac { \partial u^2 } { \partial \alpha_2} \
\frac { \partial u^2 } { \partial \alpha_3} \
u^3 \
\frac { \partial u^3 } { \partial \alpha_1} \
\frac { \partial u^3 } { \partial \alpha_2} \
\frac { \partial u^3 } { \partial \alpha_3} \
\end{array} 
\right) 
$
End of explanation
"""


E=zeros(6,9)
E[0,0]=1
E[1,4]=1
E[2,8]=1
E[3,1]=1
E[3,3]=1
E[4,2]=1
E[4,6]=1
E[5,5]=1
E[5,7]=1
E


Q=E*B
Q=simplify(Q)
Q

"""
Explanation: Deformations tensor
End of explanation
"""


T=zeros(12,6)
T[0,0]=1
T[0,2]=alpha3
T[1,1]=1
T[1,3]=alpha3
T[3,2]=1

T[8,4]=1
T[9,5]=1
T

Q=E*B*T
Q=simplify(Q)
Q

"""
Explanation: Tymoshenko theory
$u^1 \left( \alpha_1, \alpha_2, \alpha_3 \right)=u\left( \alpha_1 \right)+\alpha_3\gamma \left( \alpha_1 \right) $
$u^2 \left( \alpha_1, \alpha_2, \alpha_3 \right)=0 $
$u^3 \left( \alpha_1, \alpha_2, \alpha_3 \right)=w\left( \alpha_1 \right) $
$ \left( 
\begin{array}{c} 
u^1 \
\frac { \partial u^1 } { \partial \alpha_1} \
\frac { \partial u^1 } { \partial \alpha_2} \
\frac { \partial u^1 } { \partial \alpha_3} \
u^2 \
\frac { \partial u^2 } { \partial \alpha_1} \
\frac { \partial u^2 } { \partial \alpha_2} \
\frac { \partial u^2 } { \partial \alpha_3} \
u^3 \
\frac { \partial u^3 } { \partial \alpha_1} \
\frac { \partial u^3 } { \partial \alpha_2} \
\frac { \partial u^3 } { \partial \alpha_3} \
\end{array} 
\right) = T \cdot 
\left( 
\begin{array}{c} 
u \
\frac { \partial u } { \partial \alpha_1} \
\gamma \
\frac { \partial \gamma } { \partial \alpha_1} \
w \
\frac { \partial w } { \partial \alpha_1} \
\end{array} 
\right) $
End of explanation
"""


from sympy import MutableDenseNDimArray
C_x = MutableDenseNDimArray.zeros(3, 3, 3, 3)

for i in range(3):
    for j in range(3):        
        for k in range(3):
            for l in range(3):
                elem_index = 'C^{{{}{}{}{}}}'.format(i+1, j+1, k+1, l+1)
                el = Symbol(elem_index)
                C_x[i,j,k,l] = el
                
C_x

"""
Explanation: Elasticity tensor(stiffness tensor)
General form
End of explanation
"""


C_x_symmetry = MutableDenseNDimArray.zeros(3, 3, 3, 3)

def getCIndecies(index):
    if (index == 0):
        return 0, 0
    elif (index == 1):
        return 1, 1
    elif (index == 2):
        return 2, 2
    elif (index == 3):
        return 0, 1
    elif (index == 4):
        return 0, 2
    elif (index == 5):
        return 1, 2
    
for s in range(6):
    for t in range(s, 6):
        i,j = getCIndecies(s)
        k,l = getCIndecies(t)
        elem_index = 'C^{{{}{}{}{}}}'.format(i+1, j+1, k+1, l+1)
        el = Symbol(elem_index)
        C_x_symmetry[i,j,k,l] = el
        C_x_symmetry[i,j,l,k] = el
        C_x_symmetry[j,i,k,l] = el
        C_x_symmetry[j,i,l,k] = el
        C_x_symmetry[k,l,i,j] = el
        C_x_symmetry[k,l,j,i] = el
        C_x_symmetry[l,k,i,j] = el
        C_x_symmetry[l,k,j,i] = el

                
C_x_symmetry

"""
Explanation: Include symmetry
End of explanation
"""


C_isotropic = MutableDenseNDimArray.zeros(3, 3, 3, 3)

C_isotropic_matrix = zeros(6)

mu = Symbol('mu')
la = Symbol('lambda')

for s in range(6):
    for t in range(s, 6):
        if (s < 3 and t < 3):
            if(t != s):
                C_isotropic_matrix[s,t] = la
                C_isotropic_matrix[t,s] = la
            else:
                C_isotropic_matrix[s,t] = 2*mu+la
                C_isotropic_matrix[t,s] = 2*mu+la
        elif (s == t):
            C_isotropic_matrix[s,t] = mu
            C_isotropic_matrix[t,s] = mu
            
for s in range(6):
    for t in range(s, 6):
        i,j = getCIndecies(s)
        k,l = getCIndecies(t)
        el = C_isotropic_matrix[s, t]
        C_isotropic[i,j,k,l] = el
        C_isotropic[i,j,l,k] = el
        C_isotropic[j,i,k,l] = el
        C_isotropic[j,i,l,k] = el
        C_isotropic[k,l,i,j] = el
        C_isotropic[k,l,j,i] = el
        C_isotropic[l,k,i,j] = el
        C_isotropic[l,k,j,i] = el

                
C_isotropic

def getCalpha(C, A, q, p, s, t):
    res = S(0)
    for i in range(3):
        for j in range(3):        
            for k in range(3):
                for l in range(3):
                    res += C[i,j,k,l]*A[q,i]*A[p,j]*A[s,k]*A[t,l]
    return simplify(res)
                    


C_isotropic_alpha = MutableDenseNDimArray.zeros(3, 3, 3, 3)

for i in range(3):
    for j in range(3):        
        for k in range(3):
            for l in range(3):
                c = getCalpha(C_isotropic, A_inv, i, j, k, l)
                C_isotropic_alpha[i,j,k,l] = c


C_isotropic_matrix_alpha = zeros(6)

for s in range(6):
    for t in range(6):
        i,j = getCIndecies(s)
        k,l = getCIndecies(t)
        C_isotropic_matrix_alpha[s,t] = simplify(C_isotropic_alpha[i,j,k,l])
        
C_isotropic_matrix_alpha

"""
Explanation: Isotropic material
End of explanation
"""


C_orthotropic = MutableDenseNDimArray.zeros(3, 3, 3, 3)

C_orthotropic_matrix = zeros(6)

for s in range(6):
    for t in range(s, 6):
        elem_index = 'C^{{{}{}}}'.format(s+1, t+1)
        el = Symbol(elem_index)
        if ((s < 3 and t < 3) or t == s):
            C_orthotropic_matrix[s,t] = el
            C_orthotropic_matrix[t,s] = el
            
for s in range(6):
    for t in range(s, 6):
        i,j = getCIndecies(s)
        k,l = getCIndecies(t)
        el = C_orthotropic_matrix[s, t]
        C_orthotropic[i,j,k,l] = el
        C_orthotropic[i,j,l,k] = el
        C_orthotropic[j,i,k,l] = el
        C_orthotropic[j,i,l,k] = el
        C_orthotropic[k,l,i,j] = el
        C_orthotropic[k,l,j,i] = el
        C_orthotropic[l,k,i,j] = el
        C_orthotropic[l,k,j,i] = el

                
C_orthotropic

"""
Explanation: Orthotropic material
End of explanation
"""


def getCalpha(C, A, q, p, s, t):
    res = S(0)
    for i in range(3):
        for j in range(3):        
            for k in range(3):
                for l in range(3):
                    res += C[i,j,k,l]*A[q,i]*A[p,j]*A[s,k]*A[t,l]
    return simplify(res)
                    


C_orthotropic_alpha = MutableDenseNDimArray.zeros(3, 3, 3, 3)

for i in range(3):
    for j in range(3):        
        for k in range(3):
            for l in range(3):
                c = getCalpha(C_orthotropic, A_inv, i, j, k, l)
                C_orthotropic_alpha[i,j,k,l] = c


C_orthotropic_matrix_alpha = zeros(6)

for s in range(6):
    for t in range(6):
        i,j = getCIndecies(s)
        k,l = getCIndecies(t)
        C_orthotropic_matrix_alpha[s,t] = C_orthotropic_alpha[i,j,k,l]
        
C_orthotropic_matrix_alpha

"""
Explanation: Orthotropic material in shell coordinates
End of explanation
"""


square_int=integrate(integrate(1+alpha3/R, (alpha3, h1, h2)), (alpha1, 0, theta*R))
expand(simplify(square_int))

"""
Explanation: Square of segment
${\displaystyle A=\int_{0}^{L}\int_{h_1}^{h_2} \left( 1+\frac{\alpha_3}{R} \right) d \alpha_1 d \alpha_3}, L=R \theta$
End of explanation
"""


