Matrix generation

Init symbols for sympy



In [1]:

    
from sympy import *
from geom_util import *
from sympy.vector import CoordSys3D
N = CoordSys3D('N')
alpha1, alpha2, alpha3 = symbols("alpha_1 alpha_2 alpha_3", real = True, positive=True)
init_printing()

%matplotlib inline

%reload_ext autoreload
%autoreload 2
%aimport geom_util

Lame params



In [8]:

    
A,K = symbols('A K')
dA1 = Symbol('A_{{,1}}')
dK1 = Symbol('K_{{,1}}')
H1=A*(1+alpha3*K)
H2=S(1)
H3=S(1)

H=[H1, H2, H3]
DIM=3
dH = zeros(DIM,DIM)
dH[0,0]=dA1+alpha3*(dA1*K+A*dK1)
dH[0,2]=A*K

dH









    Out[8]:





$$\left[\begin{matrix}A_{{,1}} + \alpha_{3} \left(A K_{{,1}} + A_{{,1}} K\right) & 0 & A K\\0 & 0 & 0\\0 & 0 & 0\end{matrix}\right]$$

Metric tensor

${\displaystyle \hat{G}=\sum_{i,j} g^{ij}\vec{R}_i\vec{R}_j}$



In [9]:

    
G_up = getMetricTensorUpLame(H1, H2, H3)

${\displaystyle \hat{G}=\sum_{i,j} g_{ij}\vec{R}^i\vec{R}^j}$



In [10]:

    
G_down = getMetricTensorDownLame(H1, H2, H3)

Christoffel symbols



In [11]:

    
DIM=3

G_down_diff = MutableDenseNDimArray.zeros(DIM, DIM, DIM)
for i in range(DIM):
    for j in range(DIM):
        for k in range(DIM):
                
            G_down_diff[i,i,k]=2*H[i]*dH[i,k]
                

GK = getChristoffelSymbols2(G_up, G_down_diff, (alpha1, alpha2, alpha3))
GK









    Out[11]:





$$\left[\begin{matrix}\left[\begin{matrix}\frac{A_{{,1}} + \alpha_{3} \left(A K_{{,1}} + A_{{,1}} K\right)}{A \left(K \alpha_{3} + 1\right)} & 0 & - A^{2} K \left(K \alpha_{3} + 1\right)\\0 & 0 & 0\\\frac{K}{K \alpha_{3} + 1} & 0 & 0\end{matrix}\right] & \left[\begin{matrix}0 & 0 & 0\\0 & 0 & 0\\0 & 0 & 0\end{matrix}\right] & \left[\begin{matrix}\frac{K}{K \alpha_{3} + 1} & 0 & 0\\0 & 0 & 0\\0 & 0 & 0\end{matrix}\right]\end{matrix}\right]$$

Gradient of vector

$ \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) = B \cdot D \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) $



In [12]:

    
def row_index_to_i_j_grad(i_row):
    return i_row // 3, i_row % 3
        

B = zeros(9, 12)
B[0,1] = S(1)
B[1,2] = S(1)

B[2,3] = S(1)

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

B[6,9] = S(1)
B[7,10] = S(1)
B[8,11] = S(1)

for row_index in range(9):
    i,j=row_index_to_i_j_grad(row_index)
    B[row_index, 0] = -GK[i,j,0]
    B[row_index, 4] = -GK[i,j,1]
    B[row_index, 8] = -GK[i,j,2]

B









    Out[12]:





$$\left[\begin{array}{cccccccccccc}- \frac{A_{{,1}} + \alpha_{3} \left(A K_{{,1}} + A_{{,1}} K\right)}{A \left(K \alpha_{3} + 1\right)} & 1 & 0 & 0 & 0 & 0 & 0 & 0 & A^{2} K \left(K \alpha_{3} + 1\right) & 0 & 0 & 0\\0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\- \frac{K}{K \alpha_{3} + 1} & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\- \frac{K}{K \alpha_{3} + 1} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\end{array}\right]$$

Physical coordinates

$u_i=u_{[i]} H_i$



In [13]:

    
P=zeros(12,12)
P[0,0]=H[0]
P[1,0]=dH[0,0]
P[1,1]=H[0]
P[2,0]=dH[0,1]
P[2,2]=H[0]
P[3,0]=dH[0,2]
P[3,3]=H[0]

P[4,4]=H[1]
P[5,4]=dH[1,0]
P[5,5]=H[1]
P[6,4]=dH[1,1]
P[6,6]=H[1]
P[7,4]=dH[1,2]
P[7,7]=H[1]

P[8,8]=H[2]
P[9,8]=dH[2,0]
P[9,9]=H[2]
P[10,8]=dH[2,1]
P[10,10]=H[2]
P[11,8]=dH[2,2]
P[11,11]=H[2]
P=simplify(P)
P









    Out[13]:





$$\left[\begin{array}{cccccccccccc}A \left(K \alpha_{3} + 1\right) & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\A_{{,1}} + \alpha_{3} \left(A K_{{,1}} + A_{{,1}} K\right) & A \left(K \alpha_{3} + 1\right) & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & A \left(K \alpha_{3} + 1\right) & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\A K & 0 & 0 & A \left(K \alpha_{3} + 1\right) & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\end{array}\right]$$



In [14]:

    
B_P = zeros(9,9)

for i in range(3):
    for j in range(3):
        
        row_index = i*3+j
        
        B_P[row_index, row_index] = 1/(H[i]*H[j])
        


Grad_U_P = simplify(B_P*B*P)
Grad_U_P









    Out[14]:





$$\left[\begin{array}{cccccccccccc}0 & \frac{1}{A \left(K \alpha_{3} + 1\right)} & 0 & 0 & 0 & 0 & 0 & 0 & \frac{K}{K \alpha_{3} + 1} & 0 & 0 & 0\\0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & \frac{1}{A \left(K \alpha_{3} + 1\right)} & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\- \frac{K}{K \alpha_{3} + 1} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{A \left(K \alpha_{3} + 1\right)} & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\end{array}\right]$$

Strain tensor

$ \left( \begin{array}{c} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ 2\varepsilon_{12} \\ 2\varepsilon_{13} \\ 2\varepsilon_{23} \\ \end{array}

\right)

\left(E + E_{NL} \left( \nabla \vec{u} \right) \right) \cdot \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)$



In [15]:

    
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









    Out[15]:





$$\left[\begin{matrix}1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\\0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0\end{matrix}\right]$$



In [16]:

    
StrainL=simplify(E*Grad_U_P)
StrainL









    Out[16]:





$$\left[\begin{array}{cccccccccccc}0 & \frac{1}{A \left(K \alpha_{3} + 1\right)} & 0 & 0 & 0 & 0 & 0 & 0 & \frac{K}{K \alpha_{3} + 1} & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\\0 & 0 & 1 & 0 & 0 & \frac{1}{A \left(K \alpha_{3} + 1\right)} & 0 & 0 & 0 & 0 & 0 & 0\\- \frac{K}{K \alpha_{3} + 1} & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & \frac{1}{A \left(K \alpha_{3} + 1\right)} & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 1 & 0\end{array}\right]$$



In [17]:

    
def E_NonLinear(grad_u):
    N = 3

    du = zeros(N, N)

    #    print("===Deformations===")

    for i in range(N):
        for j in range(N):
            index = i*N+j
            du[j,i] = grad_u[index]

    #    print("========")
    
    I = eye(3)

    a_values = S(1)/S(2) * du * G_up


    E_NL = zeros(6,9)
    E_NL[0,0] = a_values[0,0]
    E_NL[0,3] = a_values[0,1]
    E_NL[0,6] = a_values[0,2]

    E_NL[1,1] = a_values[1,0]
    E_NL[1,4] = a_values[1,1]
    E_NL[1,7] = a_values[1,2]

    E_NL[2,2] = a_values[2,0]
    E_NL[2,5] = a_values[2,1]
    E_NL[2,8] = a_values[2,2]

    E_NL[3,1] = 2*a_values[0,0]
    E_NL[3,4] = 2*a_values[0,1]
    E_NL[3,7] = 2*a_values[0,2]

    E_NL[4,0] = 2*a_values[2,0]
    E_NL[4,3] = 2*a_values[2,1]
    E_NL[4,6] = 2*a_values[2,2]

    E_NL[5,2] = 2*a_values[1,0]
    E_NL[5,5] = 2*a_values[1,1]
    E_NL[5,8] = 2*a_values[1,2]


    return E_NL


%aimport geom_util
u=getUHat3DPlane(alpha1, alpha2, alpha3)


# u=getUHatU3Main(alpha1, alpha2, alpha3)

gradu=B*u


E_NL = E_NonLinear(gradu)*B

E_NL









    Out[17]:





$$\left[\begin{array}{cccccccccccc}- \frac{K}{K \alpha_{3} + 1} \left(- \frac{K u_{1}}{2 \left(K \alpha_{3} + 1\right)} + \frac{u_{3,1}}{2}\right) - \frac{1}{A^{3} \left(K \alpha_{3} + 1\right)^{3}} \left(A_{{,1}} + \alpha_{3} \left(A K_{{,1}} + A_{{,1}} K\right)\right) \left(\frac{K u_{3}}{2} A^{2} \left(K \alpha_{3} + 1\right) + \frac{u_{1,1}}{2} - \frac{u_{1} \left(A_{{,1}} + \alpha_{3} \left(A K_{{,1}} + A_{{,1}} K\right)\right)}{2 A \left(K \alpha_{3} + 1\right)}\right) & \frac{1}{A^{2} \left(K \alpha_{3} + 1\right)^{2}} \left(\frac{K u_{3}}{2} A^{2} \left(K \alpha_{3} + 1\right) + \frac{u_{1,1}}{2} - \frac{u_{1} \left(A_{{,1}} + \alpha_{3} \left(A K_{{,1}} + A_{{,1}} K\right)\right)}{2 A \left(K \alpha_{3} + 1\right)}\right) & 0 & 0 & 0 & 0 & 0 & 0 & \frac{K}{K \alpha_{3} + 1} \left(\frac{K u_{3}}{2} A^{2} \left(K \alpha_{3} + 1\right) + \frac{u_{1,1}}{2} - \frac{u_{1} \left(A_{{,1}} + \alpha_{3} \left(A K_{{,1}} + A_{{,1}} K\right)\right)}{2 A \left(K \alpha_{3} + 1\right)}\right) & - \frac{K u_{1}}{2 \left(K \alpha_{3} + 1\right)} + \frac{u_{3,1}}{2} & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\- \frac{K \left(- \frac{K u_{1}}{2 \left(K \alpha_{3} + 1\right)} + \frac{u_{1,3}}{2}\right)}{A^{2} \left(K \alpha_{3} + 1\right)^{3}} & 0 & 0 & \frac{- \frac{K u_{1}}{2 \left(K \alpha_{3} + 1\right)} + \frac{u_{1,3}}{2}}{A^{2} \left(K \alpha_{3} + 1\right)^{2}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{u_{3,3}}{2}\\0 & 0 & \frac{2}{A^{2} \left(K \alpha_{3} + 1\right)^{2}} \left(\frac{K u_{3}}{2} A^{2} \left(K \alpha_{3} + 1\right) + \frac{u_{1,1}}{2} - \frac{u_{1} \left(A_{{,1}} + \alpha_{3} \left(A K_{{,1}} + A_{{,1}} K\right)\right)}{2 A \left(K \alpha_{3} + 1\right)}\right) & 0 & 0 & 0 & 0 & 0 & 0 & 0 & - \frac{K u_{1}}{K \alpha_{3} + 1} + u_{3,1} & 0\\- \frac{K u_{3,3}}{K \alpha_{3} + 1} - \frac{2}{A^{3} \left(K \alpha_{3} + 1\right)^{3}} \left(A_{{,1}} + \alpha_{3} \left(A K_{{,1}} + A_{{,1}} K\right)\right) \left(- \frac{K u_{1}}{2 \left(K \alpha_{3} + 1\right)} + \frac{u_{1,3}}{2}\right) & \frac{2 \left(- \frac{K u_{1}}{2 \left(K \alpha_{3} + 1\right)} + \frac{u_{1,3}}{2}\right)}{A^{2} \left(K \alpha_{3} + 1\right)^{2}} & 0 & 0 & 0 & 0 & 0 & 0 & \frac{2 K}{K \alpha_{3} + 1} \left(- \frac{K u_{1}}{2 \left(K \alpha_{3} + 1\right)} + \frac{u_{1,3}}{2}\right) & u_{3,3} & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\end{array}\right]$$

Linear Stiffness matrix



In [18]:

    
%aimport geom_util
C_tensor = getIsotropicStiffnessTensor()
C = convertStiffnessTensorToMatrix(C_tensor)
C









    Out[18]:





$$\left[\begin{matrix}\lambda + 2 \mu & \lambda & \lambda & 0 & 0 & 0\\\lambda & \lambda + 2 \mu & \lambda & 0 & 0 & 0\\\lambda & \lambda & \lambda + 2 \mu & 0 & 0 & 0\\0 & 0 & 0 & \mu & 0 & 0\\0 & 0 & 0 & 0 & \mu & 0\\0 & 0 & 0 & 0 & 0 & \mu\end{matrix}\right]$$



In [22]:

    
StrainL.T*C*StrainL*H1*H2*H3









    Out[22]:





$$\left[\begin{array}{cccccccccccc}\frac{A K^{2} \mu}{K \alpha_{3} + 1} & 0 & 0 & - A K \mu & 0 & 0 & 0 & 0 & 0 & - \frac{K \mu}{K \alpha_{3} + 1} & 0 & 0\\0 & \frac{\lambda + 2 \mu}{A \left(K \alpha_{3} + 1\right)} & 0 & 0 & 0 & 0 & \lambda & 0 & \frac{K \left(\lambda + 2 \mu\right)}{K \alpha_{3} + 1} & 0 & 0 & \lambda\\0 & 0 & A \mu \left(K \alpha_{3} + 1\right) & 0 & 0 & \mu & 0 & 0 & 0 & 0 & 0 & 0\\- A K \mu & 0 & 0 & A \mu \left(K \alpha_{3} + 1\right) & 0 & 0 & 0 & 0 & 0 & \mu & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & \mu & 0 & 0 & \frac{\mu}{A \left(K \alpha_{3} + 1\right)} & 0 & 0 & 0 & 0 & 0 & 0\\0 & \lambda & 0 & 0 & 0 & 0 & A \left(\lambda + 2 \mu\right) \left(K \alpha_{3} + 1\right) & 0 & A K \lambda & 0 & 0 & A \lambda \left(K \alpha_{3} + 1\right)\\0 & 0 & 0 & 0 & 0 & 0 & 0 & A \mu \left(K \alpha_{3} + 1\right) & 0 & 0 & A \mu \left(K \alpha_{3} + 1\right) & 0\\0 & \frac{K \left(\lambda + 2 \mu\right)}{K \alpha_{3} + 1} & 0 & 0 & 0 & 0 & A K \lambda & 0 & \frac{A K^{2} \left(\lambda + 2 \mu\right)}{K \alpha_{3} + 1} & 0 & 0 & A K \lambda\\- \frac{K \mu}{K \alpha_{3} + 1} & 0 & 0 & \mu & 0 & 0 & 0 & 0 & 0 & \frac{\mu}{A \left(K \alpha_{3} + 1\right)} & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & A \mu \left(K \alpha_{3} + 1\right) & 0 & 0 & A \mu \left(K \alpha_{3} + 1\right) & 0\\0 & \lambda & 0 & 0 & 0 & 0 & A \lambda \left(K \alpha_{3} + 1\right) & 0 & A K \lambda & 0 & 0 & A \left(\lambda + 2 \mu\right) \left(K \alpha_{3} + 1\right)\end{array}\right]$$

Mass matrix



In [20]:

    
rho=Symbol('rho')
B_h=zeros(3,12)
B_h[0,0]=1
B_h[1,4]=1
B_h[2,8]=1
M=simplify(rho*P.T*B_h.T*G_up*B_h*P)
M









    Out[20]:





$$\left[\begin{array}{cccccccccccc}\rho & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & \rho & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \rho & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\end{array}\right]$$