# Overview

I made a 3D skeleton structure analysis program in Python that I've always wanted. In the plane frame analysis program created earlier, the rigidity matrix and coordinate transformation matrix are mainly changed to support beam elements with 1 node and 6 degrees of freedom.

~~ The test is not enough yet, but ~~ I was inspired by the article below and posted it. http://qiita.com/sasaco/items/917429a497c6e9a08401

~~ Since the program is created on Jupyter, this article is also described according to the description of my Jupyter. The analysis program by Python is divided into three parts. ~~ I'm always worried, but the file export part is too cluttered. Can't you write better? .. ..

~~ The calculation is going around. .. .. ~~

The programming environment is Machine: MacBook Pro (Retina, 13-inch, Mid 2014) OS: macOS Sierra Python version: Python 3.6.1

# Theory of three-dimensional skeleton structure analysis

## Basic matters

The theory development and program include the following conditions.

• Assuming minute deformation of linear elastic material
• Consider inertial force as physical force
• Consider strain due to temperature change as initial strain

In this program, assuming that the model will be large-scale, sparse matrix processing is introduced when solving simultaneous linear equations.

## Element stiffness equation

The general form of the element stiffness equation is the same as in the two-dimensional stress analysis.

### Equation for finding displacement

$$[\boldsymbol{k}]\{\boldsymbol{u_{nd}}\}=\{\boldsymbol{f}\}+\{\boldsymbol{f_t}\}+\{\boldsymbol{f_b}\}$$

\begin{align}
&[\boldsymbol{k}]=\int_V[\boldsymbol{B}]^T[\boldsymbol{D_e}][\boldsymbol{B}]dV  & & \text{(Element stiffness matrix)}\\
&\{\boldsymbol{f}\}=\int_A [\boldsymbol{N}]^T\{\boldsymbol{S}\}dA  & & \text{(Nodal external force vector)}\\
&\{\boldsymbol{f_t}\}=\int_V[\boldsymbol{B}]^T[\boldsymbol{D_e}]\{\boldsymbol{\epsilon_0}\}dV  & & \text{(Initial node distortion vector)}\\
&\{\boldsymbol{f_b}\}=\gamma\left(\int_V [\boldsymbol{N}]^T [\boldsymbol{N}]dV\right)\{\boldsymbol{w_{nd}}\} & & \text{(Nodal inertial force vector)}
\end{align}


### Element internal stress calculation formula

$$\{\boldsymbol{\sigma}\}=[\boldsymbol{D_e}]\{\boldsymbol{\epsilon}-\boldsymbol{\epsilon_0}\}$$

\begin{align}
&\{\boldsymbol{\epsilon}\}  & &\text{: Strain vector of the element obtained from the displacement obtained by solving the stiffness equation}\\
&\{\boldsymbol{\epsilon_0}\} & &\text{: Initial strain such as temperature strain}\\
&[\boldsymbol{D_e}]          & &\text{: Stress-strain relation matrix}
\end{align}


## Element stiffness matrix of three-dimensional skeleton structure

The element stiffness matrix of the three-dimensional skeleton structure is a 12 x 12 square matrix. This can be derived by extending the one for the plane skeleton, except to consider the twist. The derivation process is omitted, and only the results are described here.

### Element stiffness equation (simplest expression)

$$\boldsymbol{\{f\}}=\boldsymbol{[k]}\boldsymbol{\{u\}}$$


### Nodal displacement vector

$$\boldsymbol{\{u\}}= \begin{Bmatrix} u_i & v_i & w_i & \theta_{xi} & \theta_{yi} & \theta_{zi} & u_j & v_j & w_j & \theta_{xj} & \theta_{yj} & \theta_{zj} \end{Bmatrix}^T$$

\begin{align}
&u & &\text{: X-axis displacement} & &\theta_x& &\text{: X-axis rotation amount (torsion angle)}\\
&v & &\text{: Y-axis displacement} & &\theta_y& &\text{: Y-axis rotation amount}\\
&w & &\text{: Z-axis displacement} & &\theta_z& &\text{: Z-axis rotation amount}
\end{align}


### Nodal force vector

$$\boldsymbol{\{f\}}= \begin{Bmatrix} N_i & S_{yi} & S_{zi} & M_{xi} & M_{yi} & M_{zi} & N_j & S_{yj} & S_{zj} & M_{xj} & M_{yj} & M_{zj} \end{Bmatrix}^T$$

\begin{align}
&N  & &\text{: Axial force (x-axis directional force)} & &M_x& &\text{: X-axis moment (torsion moment)}\\
&S_y& &\text{: Y-axis shear force}   & &M_y& &\text{: Y-axis bending moment}\\
&S_z& &\text{: Z-axis shear force}   & &M_z& &\text{: Z-axis bending moment}
\end{align}


### Element stiffness matrix

$$\boldsymbol{[k]}= \begin{bmatrix} \cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 & -\cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 \\ 0 & \cfrac{12EI_z}{L^3} & 0 & 0 & 0 & \cfrac{6EI_z}{L^2} & 0 & -\cfrac{12EI_z}{L^3} & 0 & 0 & 0 & \cfrac{6EI_z}{L^2} \\ 0 & 0 & \cfrac{12EI_y}{L^3} & 0 & -\cfrac{6EI_y}{L^2} & 0 & 0 & 0 & -\cfrac{12EI_y}{L^3} & 0 & -\cfrac{6EI_y}{L^2} & 0 \\ 0 & 0 & 0 & \cfrac{GJ}{L} & 0 & 0 & 0 & 0 & 0 & -\cfrac{GJ}{L} & 0 & 0 \\ 0 & 0 & -\cfrac{6EI_y}{L^2} & 0 & \cfrac{4EI_y}{L} & 0 & 0 & 0 & \cfrac{6EI_y}{L^2} & 0 & \cfrac{2EI_y}{L} & 0 \\ 0 & \cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{4EI_z}{L} & 0 & -\cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{2EI_z}{L} \\ -\cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 & \cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 \\ 0 & -\cfrac{12EI_z}{L^3} & 0 & 0 & 0 & -\cfrac{6EI_z}{L^2} & 0 & \cfrac{12EI_z}{L^3} & 0 & 0 & 0 & -\cfrac{6EI_z}{L^2} \\ 0 & 0 & -\cfrac{12EI_y}{L^3} & 0 & \cfrac{6EI_y}{L^2} & 0 & 0 & 0 & \cfrac{12EI_y}{L^3} & 0 & \cfrac{6EI_y}{L^2} & 0 \\ 0 & 0 & 0 & -\cfrac{GJ}{L} & 0 & 0 & 0 & 0 & 0 & \cfrac{GJ}{L} & 0 & 0 \\ 0 & 0 & -\cfrac{6EI_y}{L^2} & 0 & \cfrac{2EI_y}{L} & 0 & 0 & 0 & \cfrac{6EI_y}{L^2} & 0 & \cfrac{4EI_y}{L} & 0 \\ 0 & \cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{2EI_z}{L} & 0 & -\cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{4EI_z}{L} \\ \end{bmatrix}$$

\begin{align}
&L& &\text{: Element length}        & &J& &\text{: Torsion constant}\\
&E& &\text{: Elastic modulus}      & &I_y& &\text{: Second moment of inertia of area around y-axis}\\
&G& &\text{: Shear modulus}& &I_z& &\text{: Z-axis geometrical moment of inertia}\\
&A& &\text{: Cross-sectional area}        & &
\end{align}


## Nodal inertial force vector of plane skeleton element

Consider the inertial force as the physical force, and give "acceleration" as the node vector. Here, the simplest method is to simply distribute the inertial force acting on the element to the nodal force. The total inertial force acting on the element is as follows as the value in the whole coordinate system.

\begin{align}
F_{m(x)}=&\frac{\gamma\cdot A \cdot L}{g}\cdot (k_x \cdot g) \\
F_{m(y)}=&\frac{\gamma\cdot A \cdot L}{g}\cdot (k_y \cdot g) \\
F_{m(z)}=&\frac{\gamma\cdot A \cdot L}{g}\cdot (k_z \cdot g) \\
\end{align}


Therefore, the nodal force $\ {\ boldsymbol {f_ {b}} }$ due to the inertial force can be determined as follows.

\begin{equation*}
\{\boldsymbol{F_{b}}\}=\{\boldsymbol{f_{b}}\}=\frac{\gamma\cdot A \cdot L}{2}
\begin{Bmatrix}
k_x \\
k_y \\
k_z \\
0   \\
0   \\
0   \\
k_x \\
k_y \\
k_z \\
0   \\
0   \\
0
\end{Bmatrix}
\end{equation*}


Here, $\ gamma$ is the unit volume weight of the element, $A$ is the element cross-sectional area, $L$ is the element length, and $g$ is the gravitational acceleration.

## Nodal temperature load vector of plane skeleton elements

It is assumed that the temperature change affects only the axial force of the member. Assuming that the temperature rise is a positive temperature change, the nodal force vector in the member coordinate system due to the temperature change is as follows.

\begin{equation*}
\{\boldsymbol{f_{t}}\}=EA\cdot\alpha\cdot\Delta T
\begin{Bmatrix}
-1 \\
0 \\
0 \\
0 \\
0 \\
0 \\
1 \\
0 \\
0 \\
0 \\
0 \\
0
\end{Bmatrix}
\end{equation*}


Here, $EA$ is the axial rigidity of the member, $\ alpha$ is the coefficient of linear expansion of the member, and $\ Delta T$ is the amount of temperature change of the member.

By transforming this coordinates, the nodal force vector in the overall coordinate system is determined.

\begin{equation*}
\{\boldsymbol{F_t}\}=[\boldsymbol{T}]^T\{\boldsymbol{f_t}\}
\end{equation*}


Here, $[\ boldsymbol {T}] ^ T$ is the transposed matrix of the coordinate transformation matrix $[\ boldsymbol {T}]$ (in this case, equal to the inverse matrix).

## Nodal external force vector of plane skeleton element

Regarding the nodal external force vector, in the program introduced here, the equivalent nodal external force is directly input from the input data file, and no special calculation processing is performed. The load is the horizontal force, vertical force, and moment in the overall coordinate system.

## Coordinate transformation matrix

The element stiffness equations obtained so far are in the "element coordinate system" fixed to the element. In an actual structure, the elements are oriented in various directions, so it is necessary to find the relationship with the unified "whole coordinate system".

Here, if the coordinate transformation matrix that converts the total coordinate system displacement to the element coordinate system displacement is $[\ boldsymbol {T}]$, The stiffness matrix $[\ boldsymbol {K}]$ in the global coordinate system and the stiffness matrix $[\ boldsymbol {k}]$ in the element coordinate system have the following relationship.

$$[\boldsymbol{K}]=[\boldsymbol{T}]^T[\boldsymbol{k}][\boldsymbol{T}]$$


This relationship can be understood from the relationship between the element coordinate system displacement $\ {\ boldsymbol {u} }$ and the global coordinate system displacement $\ {\ boldsymbol {U} }$ as follows. ..

$$\{\boldsymbol{u}\}=[\boldsymbol{T}]\{\boldsymbol{U}\} \qquad\qquad \{\boldsymbol{u}\}^T=\{\boldsymbol{U}\}^T [\boldsymbol{T}]^T$$


Here, $[\ boldsymbol {T}]$ is a coordinate conversion matrix that converts the displacement / nodal force of the entire coordinate system into the displacement / nodal force of the element coordinate system.

### Specific representation of the coordinate transformation matrix

#### Direction cosine of member axis (x axis)

Let the nodes that make up the element be $i$ and $j$, and let their coordinates be $(x_i, y_i, z_i)$, $(x_j, y_j, z_j)$.

\begin{gather}
\ell=\sqrt{(x_j-x_i)^2+(y_j-y_i)^2+(z_j-z_i)^2}
\end{gather}


#### Coordinate transformation submatrix

$\ boldsymbol {[t]}$ is a submatrix (3 x 3) of the coordinate transformation matrix (12 x 12) that transforms the entire coordinate system into the member coordinate system.

$$\boldsymbol{\{x\}}=\boldsymbol{[t]}\boldsymbol{\{X\}}$$

\begin{align}
&\boldsymbol{\{x\}}=\{x,y,z\}^T  & &\text{: Local coordinates}\\
&\boldsymbol{\{X\}}＝\{X,Y,Z\}^T & &\text{: Global coordinate}
\end{align}

• When the member long axis (member coordinate system x axis) is not parallel to the overall coordinate system Z axis
$$\boldsymbol{[t]}= \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos\theta & \sin\theta \\ 0 & -\sin\theta & \cos\theta \end{bmatrix} \cdot \begin{bmatrix} l & m & n \\ -\cfrac{m}{\sqrt{l^2+m^2}} & \cfrac{l}{\sqrt{l^2+m^2}} & 0 \\ -\cfrac{l \cdot n}{\sqrt{l^2+m^2}} & -\cfrac{m \cdot n}{\sqrt{l^2+m^2}} & \sqrt{l^2+m^2} \end{bmatrix} \qquad \text{(\theta : chord angle)}$$

• When the member long axis (member coordinate system x axis) is parallel to the overall coordinate system Z axis
At this time, $x_j-x_i = 0$, $y_j-y_i = 0$, and $l = 0$, $m = 0$, so it is necessary to rewrite $\ boldsymbol {[t]}$. is there.
$$\boldsymbol{[t]}= \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos\theta & \sin\theta \\ 0 & -\sin\theta & \cos\theta \end{bmatrix} \cdot \begin{bmatrix} 0 & 0 & n \\ n & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix} \qquad \text{(\theta : chord angle)}$$


#### Coordinate transformation matrix

$$\boldsymbol{[T]}= \begin{bmatrix} \boldsymbol{[t]} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} \\\ \boldsymbol{0} & \boldsymbol{[t]} & \boldsymbol{0} & \boldsymbol{0} \\\ \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{[t]} & \boldsymbol{0} \\\ \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} &\boldsymbol{[t]} \end{bmatrix}$$


## Overall stiffness equation

From the above, the stiffness equation in the overall coordinate system is expressed as follows.

$$[\boldsymbol{K}]\{\boldsymbol{U}\}=\{\boldsymbol{F}\}+\{\boldsymbol{F_b}\}+\{\boldsymbol{F_t}\}$$

\begin{align}
&\{\boldsymbol{U}\}                          & & \text{Nodal displacement vector in global coordinate system}\\
&\{\boldsymbol{F}\}                          & & \text{Node external force vector in the global coordinate system}\\
&\{\boldsymbol{F_b}\}=\{\boldsymbol{f_b}\}   & & \text{Nodal inertial force vector in global coordinate system}\\
&\{\boldsymbol{F_t}\}=[\boldsymbol{T}]^T\{\boldsymbol{f_t}\}  & & \text{Nodal temperature load vector in global coordinate system}\\
&[\boldsymbol{K}]=[\boldsymbol{T}]^T[\boldsymbol{k}][\boldsymbol{T}]   & & \text{Rigidity matrix in the global coordinate system}\\
&[\boldsymbol{k}]                            & & \text{Rigidity matrix in element coordinate system}\\
&\{\boldsymbol{f_b}\}                          & & \text{Nodal inertial force vector in element coordinate system}\\
&\{\boldsymbol{f_t}\}                          & & \text{Nodal temperature load vector in element coordinate system}\\
&[\boldsymbol{T}]                            & & \text{Coordinate transformation matrix}
\end{align}


By solving the above equation, the node displacement $\ {\ boldsymbol {U} }$ in the global coordinate system can be obtained.

## Element cross-sectional force

For the element section force, the total coordinate system displacement $\ {\ boldsymbol {U} }$ obtained above is converted to the element coordinate system displacement by the coordinate conversion matrix $[\ boldsymbol {T}]$. Calculated by multiplying the element stiffness matrix $[\ boldsymbol {k}]$. That is, the element section force $\ {\ boldsymbol {f_ {sec}} }$ is calculated by the following.

$$\{\boldsymbol{f_{sec}}\}=[\boldsymbol{k}]\cdot ([\boldsymbol{T}]\{\boldsymbol{U}\})$$


Note that the calculated axial force calculated above does not include the effect of temperature change, so it is necessary to correct the temperature effect of the axial force as follows.

\begin{align}
N_1'=&N_1+EA\cdot\alpha\cdot\Delta T \\
N_2'=&N_2-EA\cdot\alpha\cdot\Delta T
\end{align}


Here, $N_1'$ and $N_2'$ are the axial forces at nodes 1 and 2 considering the temperature change, and $N_1$ and $N_2$ are the axial forces calculated from the displacement as the solution of the simultaneous equations. $EA$ is the axial stiffness of the member, $\ alpha$ is the coefficient of linear expansion of the member, and $\ Delta T$ is the amount of temperature change of the member.

# Three-dimensional frame structure analysis program

## Program overview

• Programming language is Python
• Perform 3D skeleton elasticity analysis The + element is a beam element with 1 element, 2 nodes, 1 node and 6 degrees of freedom.
• In the global coordinate system, the Cartesian coordinate system on the horizontal plane has the x-axis positive direction in the right direction and the y-axis positive direction in the upward direction.
• Can handle node external force, node forced displacement (including zero displacement), node temperature change, and inertial force on elements as loads.
• External force can handle only those that act on nodes. Cannot handle the load in the middle of the element
• Simultaneous linear equations are solved by introducing forced displacement of nodes without changing their size (1 degree of freedom of nodes x total number of nodes).
• Simultaneous linear equations are solved by scipy.sparse.linalg.spsolve () </ var> after CSR conversion.
• Output is node displacement and element section force
• The character delimiter of the input data file and output data file is blank.
• Enter the input data file name and output data file name as command line arguments.

## Explanation of boundary conditions that can be handled

Item Description
Nodal temperature change Specify the amount of temperature change at all nodes. Let the temperature rise be the positive value
Element inertial force In the material property data, specify the acceleration in the x / y / z direction as the ratio to the gravitational acceleration
Nodal forced displacement Specify the number of nodes, nodal number, and displacement value to give forced displacement

## Program configuration (program name: py_fem_3dfrmx.py)

   1.Module loading
2.Function: data entry(def inpdata_3dfrm)
(2)Array declaration
(4)Read element composition node number and material number
(5)Read node coordinates and node temperature change
3.Function: Write input data(def prinp_3dfrm)
4.Function: Export calculation result(def prout_3dfrm)
5.Function: Create element stiffness matrix(def sm_3dfrm)
6.Function: Create coordinate conversion matrix(def tm_3dfrm)
7.Function: Node inertial force vector creation(def bfvec_3dfrm)
8.Function: Node temperature load vector creation(def tfvec_3dfrm)
9.Function: Element section force calculation(def calsecf_3dfrm)
10.Function: Main routine(def main_3dfrm)
(1)Get I / O file name from command line
(2)Data entry
(3)Input data export
(4)Rigid matrix assembly
(5)Boundary condition processing (processing of known displacement degree)
(6)Solving the stiffness equation (displacement calculation)
(7)Boundary condition processing (incorporation of known displacement)
(8)Element section force calculation
(9)Export calculation results
(10)Information display such as calculation time
11.Main routine execution


## Execution command and I / O data

### Analysis execution command

python3 py_fem_3dfrmx.py inp.txt out.txt

Item Description
py_fem_3dfrmx.py Python FEM program name
inp.txt Input data file name (blank delimited data)
out.txt Output data file name (blank separated data)

### Input data file format

npoin  nele  nsec  npfix  nlod                   # Basic values for analysis
E  po  A  Ix  Iy  Iz  theta  alpha  gamma  gkX  gkY  gkZ
..... (1 to nsec) .....                      # Material properties
node_1  node_2  isec                             # Element connectivity, material set number
..... (1 to nele) .....                      #
x  y  z  deltaT                                  # Node coordinate, temperature change of node
..... (1 to npoin) .....                     #
lp  kox  koy  koz  kmx  kmy  kmz  rdis_x  rdis_y  rdis_z  rrot_x  rrot_y  rrot_z
..... (1 to npfix) .....                     # Restricted node and restrict conditions
..... (1 to nlod) .....                      #     (omit data input if nlod=0)

Item Description
npoin, nele, nsec Number of nodes, number of elements, number of cross-sectional characteristics
npfix, nlod Number of constrained nodes, number of loaded nodes
E, po, A Element elastic modulus, element Poisson's ratio, element cross-sectional area
Ix, Iy, Iz Torsion constant, moment of inertia of area around y-axis, moment of inertia of area around z-axis >
theta, alpha Code angle, coefficient of linear expansion
gamma, gkX, gkY, gkZ Unit volume weight, acceleration in x / y / z direction (ratio of g)
node_1, node_2, isec Node 1, Node 2, Sectional characteristic number
x, y, z, deltaT Node x coordinate, node y coordinate, node z coordinate, node temperature change
lp, kox, koy, koz Constraint node number, x / y / z Directional displacement With or without constraint (constraint: 1, freedom: 0)
kmx, kmy, kmz x / y / z With or without axial rotation constraint (constraint: 1, freedom: 0) < / tr>
r_disx, r_disy, rdis_z x / y / z direction known displacement (enter 0 even if unconstrained)
rrot_x, rrot_y, rrot_z x / y / z Known rotation amount around axis (enter 0 even if unconstrained)
mp_x, mp_y, mp_z x-axis moment (torsion), y-axis moment, z-axis moment

### Output data file format

npoin  nele  nsec npfix  nlod
(Each value of above)
sec  E      po     A    J    Iy    Iz    theta
sec  alpha  gamma  gkX  gkY  gkZ
sec   : Material number
E     : Elastic modulus
po    : Poisson's ratio
A     : Section area
J     ; Torsional constant
Iy    : Moment of inertia around y-axis
Iz    : Moment of inertia around z-axis
theta : Chord angle
alpha : Thermal expansion coefficient
gamma : Unit weight
gkX   : Acceleration in X-direction (ratio to g)
gkY   : Acceleration in Y-direction (ratio to g)
gkZ   : Acceleration in Z-direction (ratio to g)
..... (1 to nsec) .....
node   x   y   z   fx   fy   fz   mx   my   mz  deltaT
node   : Node number
x      : x-coordinate
y      : y-coordinate
mx     : Moment load around X-axis
my     : Moment load around Y-axis
mz     : Moment load around Z-axis
deltaT : Temperature change of node
..... (1 to npoin) .....
node   kox   koy   koz   kmx   kmy   kmz   rdis_x   rdis_y   rdis_z   rrot_x   rrot_y   rrot_z
node    : Node number
kox     : Index of restriction in X-direction (0: free, 1: fixed)
koy     : Index of restriction in Y-direction (0: free, 1: fixed)
koz     : Index of restriction in Z-direction (0: free, 1: fixed)
kmx     : Index of restriction in rotation around X-axis  (0: free, 1: fixed)
kmy     : Index of restriction in rotation around Y-axis  (0: free, 1: fixed)
kmz     : Index of restriction in rotation around Z-axis  (0: free, 1: fixed)
rdis_x  : Given displacement in X-direction
rdis_y  : Given displacement in Y-direction
rdis_z  : Given displacement in Z-direction
rrot_x  : Given displacement in rotation around X-axis
rrot_y  : Given displacement in rotation around Y-axis
rrot_z  : Given displacement in rotation around Z-axis
..... (1 to npfix) .....
elem     i     j   sec
elem : Element number
i    : Node number of start point
j    : Node number of end point
sec  : Material number
..... (1 to nele) .....
node   dis-x   dis-y   dis-z   rot-x   rot-y   rot-z
node  : Node number
dis-x : Displacement in X-direction
dis-y : Displacement in Y-direction
dis-z : Displacement in Z-direction
rot-x : Rotation around X-axis
rot-y : Rotation around Y-axis
rot-z : Rotation around Z-axis
..... (1 to npoin) .....
elem nodei   N_i   Sy_i   Sz_i   Mx_i   My_i   Mz_i
elem nodej   N_j   Sy_j   Sz_j   Mx_j   My_j   Mz_j
elem  : Element number
nodei : First node nymber
nodej : Second node number
N_i   : Axial force of node-i
Sy_i  : Shear force of node-i in y-direction
Sz_i  : Shear force of node-i in z-direction
Mx_i  : Moment of node-i around x-axis (torsional moment)
My_i  : Moment of node-i around y-axis (bending moment)
Mz_i  : Moment of node-i around z-axis (bending moment)
N_j   : Axial force of node-j
Sy_j  : Shear force of node-j in y-direction
Sz_j  : Shear force of node-j in z-direction
Mx_j  : Moment of node-j around x-axis (torsional moment)
My_j  : Moment of node-j around y-axis (bending moment)
Mz_j  : Moment of node-j around z-axis (bending moment)
..... (1 to nele) .....
n=(total degrees of freedom)  time=(calculation time)

Item Description
dis-x, dis-y, dis-z x-direction displacement, y-direction displacement, z-direction displacement
rot-x, rot-y, rot-z x-axis rotation amount, y-axis rotation amount, z-axis rotation amount td>
N, Sy, Sz Axial force, y-direction shear force, z-direction shear force
Mx, My, Mz x-axis moment (torsion), y-axis bending moment, z-axis bending moment

## programming

The main body of the program is shown below according to the "program structure".

• numpy: For array operation
• scipy.sparse.linalg.spsolve: Solving sparse matrix simultaneous linear equations
• scipy.sparse import csr_matrix: Sparse matrix compression storage
• sys: use command line arguments
• time: For calculation time measurement
# ======================
# 3D Frame Analysis
# ======================
import numpy as np
from scipy.sparse.linalg import spsolve
from scipy.sparse import csr_matrix
import sys
import time


### 2. Function: Data entry (def inpdata_3dfrm)

• Read basic numbers ( npoin, nele, nsec, npfix, nlod </ var>)
• Array declaration ( ae, node, x, deltaT, mpfix, rdis, fp </ var>)
• Read material property data (array ae </ var>)
• Read element constituent node numbers and material numbers (array node </ var>)
• Read node coordinates and node temperature change (array x, deltaT </ var>)
• Boundary condition read (array mpfix, rdis </ var>)
def inpdata_3dfrm(fnameR,nod,nfree):
f=open(fnameR,'r')
text=text.strip()
text=text.split()
npoin=int(text[0]) # Number of nodes
nele =int(text[1]) # Number of elements
nsec =int(text[2]) # Number of sections
npfix=int(text[3]) # Number of restricted nodes
nlod =int(text[4]) # Number of loaded nodes
# array declaration
ae    =np.zeros((12,nsec),dtype=np.float64)      # Section characteristics
node  =np.zeros((nod+1,nele),dtype=np.int)       # node-element relationship
x     =np.zeros((3,npoin),dtype=np.float64)      # Coordinates of nodes
deltaT=np.zeros(npoin,dtype=np.float64)          # Temperature change of nodes
mpfix =np.zeros((nfree,npoin),dtype=np.int)      # Ristrict conditions
rdis  =np.zeros((nfree,npoin),dtype=np.float64)  # Ristricted displacement
fp    =np.zeros(nfree*npoin,dtype=np.float64)    # External force vector
# section characteristics
for i in range(0,nsec):
text=text.strip()
text=text.split()
ae[0,i] =float(text[0])  # E     : elastic modulus
ae[1,i] =float(text[1])  # po    : Poisson's ratio
ae[2,i] =float(text[2])  # aa    : section area
ae[3,i] =float(text[3])  # aix   : tortional constant
ae[4,i] =float(text[4])  # aiy   : moment of inertia aroutd y-axis
ae[5,i] =float(text[5])  # aiz   : moment of inertia around z-axis
ae[6,i] =float(text[6])  # theta : chord angle
ae[7,i] =float(text[7])  # alpha : thermal expansion coefficient
ae[8,i] =float(text[8])  # gamma : unit weight of material
ae[9,i] =float(text[9])  # gkX   : acceleration in X-direction
ae[10,i]=float(text[10]) # gkY   : acceleration in Y-direction
ae[11,i]=float(text[11]) # gkZ   : acceleration in Z-direction
# element-node
for i in range(0,nele):
text=text.strip()
text=text.split()
node[0,i]=int(text[0]) #node_1
node[1,i]=int(text[1]) #node_2
node[2,i]=int(text[2]) #section characteristic number
# node coordinates
for i in range(0,npoin):
text=text.strip()
text=text.split()
x[0,i]=float(text[0])    # x-coordinate
x[1,i]=float(text[1])    # y-coordinate
x[2,i]=float(text[2])    # z-coordinate
deltaT[i]=float(text[3]) # Temperature change
# boundary conditions (0:free, 1:restricted)
for i in range(0,npfix):
text=text.strip()
text=text.split()
lp=int(text[0])              #fixed node
mpfix[0,lp-1]=int(text[1])   #fixed in x-direction
mpfix[1,lp-1]=int(text[2])   #fixed in y-direction
mpfix[2,lp-1]=int(text[3])   #fixed in z-direction
mpfix[3,lp-1]=int(text[4])   #fixed in rotation around x-axis
mpfix[4,lp-1]=int(text[5])   #fixed in rotation around y-axis
mpfix[5,lp-1]=int(text[6])   #fixed in rotation around z-axis
rdis[0,lp-1]=float(text[7])  #fixed displacement in x-direction
rdis[1,lp-1]=float(text[8])  #fixed displacement in y-direction
rdis[2,lp-1]=float(text[9])  #fixed displacement in z-direction
rdis[3,lp-1]=float(text[10]) #fixed rotation around x-axis
rdis[4,lp-1]=float(text[11]) #fixed rotation around y-axis
rdis[5,lp-1]=float(text[12]) #fixed rotation around z-axis
if 0&lt;nlod:
for i in range(0,nlod):
text=text.strip()
text=text.split()
fp[6*lp-3]=float(text[4]) #moment around x-axis
fp[6*lp-2]=float(text[5]) #moment around y-axis
fp[6*lp-1]=float(text[6]) #moment around z-axis
f.close()
return npoin,nele,nsec,npfix,nlod,ae,node,x,deltaT,mpfix,rdis,fp


### 3. Function: Write input data (def prinp_3dfrm)

def prinp_3dfrm(fnameW,npoin,nele,nsec,npfix,nlod,ae,node,x,deltaT,mpfix,rdis,fp):
fout=open(fnameW,'w')
# print out of input data
print('{0:&gt;5s} {1:&gt;5s} {2:&gt;5s} {3:&gt;5s} {4:&gt;5s}'.format('npoin','nele','nsec','npfix','nlod'),file=fout)
print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d}'.format(npoin,nele,nsec,npfix,nlod),file=fout)
print('{0:&gt;5s} {1:&gt;15s} {2:&gt;15s} {3:&gt;15s} {4:&gt;15s} {5:&gt;15s} {6:&gt;15s} {7:&gt;15s}'
.format('sec','E','po','A','J','Iy','Iz','theta'),file=fout)
print('{0:&gt;5s} {1:&gt;15s} {2:&gt;15s} {3:&gt;15s} {4:&gt;15s} {5:&gt;15s}'
.format('sec','alpha','gamma','gkX','gkY','gkZ'),file=fout)
for i in range(0,nsec):
print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e} {7:15.7e}'
.format(i+1,ae[0,i],ae[1,i],ae[2,i],ae[3,i],ae[4,i],ae[5,i],ae[6,i]),file=fout)
print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e}'
.format(i+1,ae[7,i],ae[8,i],ae[9,i],ae[10,i],ae[11,i]),file=fout)
print('{0:&gt;5s} {1:&gt;15s} {2:&gt;15s} {3:&gt;15s} {4:&gt;15s} {5:&gt;15s} {6:&gt;15s} {7:&gt;15s} {8:&gt;15s} {9:&gt;15s} {10:&gt;15s}'
.format('node','x','y','z','fx','fy','fz','mx','my','mz','deltaT'),file=fout)
for i in range(0,npoin):
lp=i+1
print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e} {7:15.7e} {8:15.7e} {9:15.7e} {10:15.7e}'
.format(lp,x[0,i],x[1,i],x[2,i],fp[6*lp-6],fp[6*lp-5],fp[6*lp-4],fp[6*lp-3],fp[6*lp-2],fp[6*lp-1],deltaT[i]),file=fout)
print('{0:&gt;5s} {1:&gt;5s} {2:&gt;5s} {3:&gt;5s} {4:&gt;5s} {5:&gt;5s} {6:&gt;5s} {7:&gt;15s} {8:&gt;15s} {9:&gt;15s} {10:&gt;15s} {11:&gt;15s} {12:&gt;15s}'
.format('node','kox','koy','koz','kmx','kmy','kmz','rdis_x','rdis_y','rdis_z','rrot_x','rrot_y','rrot_z'),file=fout)
for i in range(0,npoin):
if 0&lt;mpfix[0,i]+mpfix[1,i]+mpfix[2,i]+mpfix[3,i]+mpfix[4,i]+mpfix[5,i]:
lp=i+1
print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d} {5:5d} {6:5d} {7:15.7e} {8:15.7e} {9:15.7e} {10:15.7e} {11:15.7e} {12:15.7e}'
.format(lp,mpfix[0,i],mpfix[1,i],mpfix[2,i],mpfix[3,i],mpfix[4,i],mpfix[5,i],rdis[0,i],rdis[1,i],rdis[2,i],rdis[3,i],rdis[4,i],rdis[5,i]),file=fout)
print('{0:&gt;5s} {1:&gt;5s} {2:&gt;5s} {3:&gt;5s}'.format('elem','i','j','sec'),file=fout)
for ne in range(0,nele):
print('{0:5d} {1:5d} {2:5d} {3:5d}'.format(ne+1,node[0,ne],node[1,ne],node[2,ne]),file=fout)
fout.close()


### 4. Function: Export calculation result (def prout_3dfrm)

def prout_3dfrm(fnameW,npoin,nele,node,disg,fsec):
fout=open(fnameW,'a')
# displacement
print('{0:&gt;5s} {1:&gt;15s} {2:&gt;15s} {3:&gt;15s} {4:&gt;15s} {5:&gt;15s} {6:&gt;15s}'
.format('node','dis-x','dis-y','dis-z','rot-x','rot-y','rot-z'),file=fout)
for i in range(0,npoin):
lp=i+1
print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e}'
.format(lp,disg[6*lp-6],disg[6*lp-5],disg[6*lp-4],disg[6*lp-3],disg[6*lp-2],disg[6*lp-1]),file=fout)
# section force
print('{0:&gt;5s} {1:&gt;5s} {2:&gt;15s} {3:&gt;15s} {4:&gt;15s} {5:&gt;15s} {6:&gt;15s} {7:&gt;15s}'
.format('elem','nodei','N_i','Sy_i','Sz_i','Mx_i','My_i','Mz_i'),file=fout)
print('{0:&gt;5s} {1:&gt;5s} {2:&gt;15s} {3:&gt;15s} {4:&gt;15s} {5:&gt;15s} {6:&gt;15s} {7:&gt;15s}'
.format('elem','nodej','N_j','Sy_j','Sz_j','Mx_j','My_j','Mz_j'),file=fout)
for ne in range(0,nele):
print('{0:5d} {1:5d} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e} {7:15.7e}'
.format(ne+1,node[0,ne],fsec[0,ne],fsec[1,ne],fsec[2,ne],fsec[3,ne],fsec[4,ne],fsec[5,ne]),file=fout)
print('{0:5d} {1:5d} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e} {7:15.7e}'
.format(ne+1,node[1,ne],fsec[6,ne],fsec[7,ne],fsec[8,ne],fsec[9,ne],fsec[10,ne],fsec[11,ne]),file=fout)
fout.close()


### 5. Function: Create element stiffness matrix (def sm_3dfrm)

$$\boldsymbol{[k]}= \begin{bmatrix} \cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 & -\cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 \\ 0 & \cfrac{12EI_z}{L^3} & 0 & 0 & 0 & \cfrac{6EI_z}{L^2} & 0 & -\cfrac{12EI_z}{L^3} & 0 & 0 & 0 & \cfrac{6EI_z}{L^2} \\ 0 & 0 & \cfrac{12EI_y}{L^3} & 0 & -\cfrac{6EI_y}{L^2} & 0 & 0 & 0 & -\cfrac{12EI_y}{L^3} & 0 & -\cfrac{6EI_y}{L^2} & 0 \\ 0 & 0 & 0 & \cfrac{GJ}{L} & 0 & 0 & 0 & 0 & 0 & -\cfrac{GJ}{L} & 0 & 0 \\ 0 & 0 & -\cfrac{6EI_y}{L^2} & 0 & \cfrac{4EI_y}{L} & 0 & 0 & 0 & \cfrac{6EI_y}{L^2} & 0 & \cfrac{2EI_y}{L} & 0 \\ 0 & \cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{4EI_z}{L} & 0 & -\cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{2EI_z}{L} \\ -\cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 & \cfrac{EA}{L} & 0 & 0 & 0 & 0 & 0 \\ 0 & -\cfrac{12EI_z}{L^3} & 0 & 0 & 0 & -\cfrac{6EI_z}{L^2} & 0 & \cfrac{12EI_z}{L^3} & 0 & 0 & 0 & -\cfrac{6EI_z}{L^2} \\ 0 & 0 & -\cfrac{12EI_y}{L^3} & 0 & \cfrac{6EI_y}{L^2} & 0 & 0 & 0 & \cfrac{12EI_y}{L^3} & 0 & \cfrac{6EI_y}{L^2} & 0 \\ 0 & 0 & 0 & -\cfrac{GJ}{L} & 0 & 0 & 0 & 0 & 0 & \cfrac{GJ}{L} & 0 & 0 \\ 0 & 0 & -\cfrac{6EI_y}{L^2} & 0 & \cfrac{2EI_y}{L} & 0 & 0 & 0 & \cfrac{6EI_y}{L^2} & 0 & \cfrac{4EI_y}{L} & 0 \\ 0 & \cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{2EI_z}{L} & 0 & -\cfrac{6EI_z}{L^2} & 0 & 0 & 0 & \cfrac{4EI_z}{L} \\ \end{bmatrix}$$

def sm_3dfrm(EA,GJ,EIy,EIz,x1,y1,z1,x2,y2,z2):
ek=np.zeros((12,12),dtype=np.float64) # local stiffness matrix
xx=x2-x1
yy=y2-y1
zz=z2-z1
el=np.sqrt(xx**2+yy**2+zz**2)
ek[ 0, 0]= EA/el
ek[ 0, 6]=-EA/el
ek[ 1, 1]= 12*EIz/el**3
ek[ 1, 5]=  6*EIz/el**2
ek[ 1, 7]=-12*EIz/el**3
ek[ 1,11]=  6*EIz/el**2
ek[ 2, 2]= 12*EIy/el**3
ek[ 2, 4]= -6*EIy/el**2
ek[ 2, 8]=-12*EIy/el**3
ek[ 2,10]= -6*EIy/el**2
ek[ 3, 3]= GJ/el
ek[ 3, 9]=-GJ/el
ek[ 4, 2]= -6*EIy/el**2
ek[ 4, 4]=  4*EIy/el
ek[ 4, 8]=  6*EIy/el**2
ek[ 4,10]=  2*EIy/el
ek[ 5, 1]=  6*EIz/el**2
ek[ 5, 5]=  4*EIz/el
ek[ 5, 7]= -6*EIz/el**2
ek[ 5,11]=  2*EIz/el
ek[ 6, 0]=-EA/el
ek[ 6, 6]= EA/el
ek[ 7, 1]=-12*EIz/el**3
ek[ 7, 5]= -6*EIz/el**2
ek[ 7, 7]= 12*EIz/el**3
ek[ 7,11]= -6*EIz/el**2
ek[ 8, 2]=-12*EIy/el**3
ek[ 8, 4]=  6*EIy/el**2
ek[ 8, 8]= 12*EIy/el**3
ek[ 8,10]=  6*EIy/el**2
ek[ 9, 3]=-GJ/el
ek[ 9, 9]= GJ/el
ek[10, 2]= -6*EIy/el**2
ek[10, 4]=  2*EIy/el
ek[10, 8]=  6*EIy/el**2
ek[10,10]=  4*EIy/el
ek[11, 1]=  6*EIz/el**2
ek[11, 5]=  2*EIz/el
ek[11, 7]= -6*EIz/el**2
ek[11,11]=  4*EIz/el
return ek


### 6. Function: Create coordinate transformation matrix (def tm_3dfrm)

\begin{gather}
\ell=\sqrt{(x_j-x_i)^2+(y_j-y_i)^2+(z_j-z_i)^2}\\
\end{gather}

$$\boldsymbol{[t_1]}= \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos\theta & \sin\theta \\ 0 & -\sin\theta & \cos\theta \end{bmatrix} \qquad \text{(\theta : chord angle)}$$

\begin{align}
&\boldsymbol{[t_2]}=
\begin{bmatrix}
0 & 0 & n \\
n & 0 & 0 \\
0 & 1 & 0
\end{bmatrix}
& &\text{When the member long axis (member coordinate system x-axis) is parallel to the overall coordinate system Z-axis}\\
&\boldsymbol{[t_2]}=
\begin{bmatrix}
l & m & n \\
-\cfrac{m}{\sqrt{l^2+m^2}} & \cfrac{l}{\sqrt{l^2+m^2}} & 0 \\
-\cfrac{l \cdot n}{\sqrt{l^2+m^2}} & -\cfrac{m \cdot n}{\sqrt{l^2+m^2}} & \sqrt{l^2+m^2}
\end{bmatrix}
& &\text{When the member long axis (member coordinate system x-axis) is not parallel to the overall coordinate system Z-axis}
\end{align}

$$\boldsymbol{[t]}=\boldsymbol{[t_1]}\cdot\boldsymbol{[t_2]}$$

$$\boldsymbol{[T]}= \begin{bmatrix} \boldsymbol{[t]} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{[t]} & \boldsymbol{0} & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{[t]} & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} &\boldsymbol{[t]} \end{bmatrix}$$

def tm_3dfrm(theta,x1,y1,z1,x2,y2,z2):
tt=np.zeros((12,12),dtype=np.float64) # transformation matrix
t1=np.zeros((3,3),dtype=np.float64)
t2=np.zeros((3,3),dtype=np.float64)
xx=x2-x1
yy=y2-y1
zz=z2-z1
el=np.sqrt(xx**2+yy**2+zz**2)
t1[0,0]=1
t1[1,1]= np.cos(theta)
t1[1,2]= np.sin(theta)
t1[2,1]=-np.sin(theta)
t1[2,2]= np.cos(theta)
ll=(x2-x1)/el
mm=(y2-y1)/el
nn=(z2-z1)/el
if x2-x1==0.0 and y2-y1==0.0:
t2[0,2]=nn
t2[1,0]=nn
t2[2,1]=1.0
else:
qq=np.sqrt(ll**2+mm**2)
t2[0,0]=ll
t2[0,1]=mm
t2[0,2]=nn
t2[1,0]=-mm/qq
t2[1,1]= ll/qq
t2[2,0]=-ll*nn/qq
t2[2,1]=-mm*nn/qq
t2[2,2]=qq
t3=np.dot(t1,t2)
tt[0:3,0:3]  =t3[0:3,0:3]
tt[3:6,3:6]  =t3[0:3,0:3]
tt[6:9,6:9]  =t3[0:3,0:3]
tt[9:12,9:12]=t3[0:3,0:3]
return tt


### 7. Function: Node inertial force vector creation (def bfvec_3dfrm)

\begin{equation*}
\{\boldsymbol{F_{b}}\}=\{\boldsymbol{f_{b}}\}=\frac{\gamma\cdot A \cdot L}{2}
\begin{Bmatrix}
k_x \\
k_y \\
k_z \\
0   \\
0   \\
0   \\
k_x \\
k_y \\
k_z \\
0   \\
0   \\
0
\end{Bmatrix}
\end{equation*}

def bfvec_3dfrm(A,gamma,gkX,gkY,gkZ,x1,y1,z1,x2,y2,z2):
# Body force vector in global coordinate system
bfe_g=np.zeros(12,dtype=np.float64)
xx=x2-x1
yy=y2-y1
zz=z2-z1
el=np.sqrt(xx**2+yy**2+zz**2)
bfe_g[0]=0.5*gamma*A*el*gkX
bfe_g[1]=0.5*gamma*A*el*gkY
bfe_g[2]=0.5*gamma*A*el*gkZ
bfe_g[6]=0.5*gamma*A*el*gkX
bfe_g[7]=0.5*gamma*A*el*gkY
bfe_g[8]=0.5*gamma*A*el*gkZ
return bfe_g


### 8. Function: Node temperature load vector creation (def tfvec_3dfrm)

\begin{equation*}
\{\boldsymbol{f_{t}}\}=EA\cdot\alpha\cdot\Delta T
\begin{Bmatrix}
-1 \\
0 \\
0 \\
0 \\
0 \\
0 \\
1 \\
0 \\
0 \\
0 \\
0 \\
0
\end{Bmatrix}
\end{equation*}

def tfvec_3dfrm(EA,alpha,tem):
# Thermal load vector  in local coordinate system
tfe_l=np.zeros(12,dtype=np.float64)
tfe_l[0]=-EA*alpha*tem
tfe_l[6]= EA*alpha*tem
return tfe_l


### 9. Function: Element section force calculation (def calsecf_3dfrm)

$N_1'$ and $N_2'$ are axial force corrections due to temperature changes.

\begin{align}
&\{\boldsymbol{f_{sec}}\}=[\boldsymbol{k}]\cdot ([\boldsymbol{T}]\{\boldsymbol{U}\})\\
&N_1'=N_1+EA\cdot\alpha\cdot\Delta T \\
&N_2'=N_2-EA\cdot\alpha\cdot\Delta T
\end{align}

def calsecf_3dfrm(EA,GJ,EIy,EIz,theta,alpha,tem,dis,x1,y1,z1,x2,y2,z2):
ek=sm_3dfrm(EA,GJ,EIy,EIz,x1,y1,z1,x2,y2,z2) # Stiffness matrix in local coordinate
tt=tm_3dfrm(theta,x1,y1,z1,x2,y2,z2) # Transformation matrix
secf=np.dot(ek,np.dot(tt,dis))
secf[0]=secf[0]+EA*alpha*tem
secf[6]=secf[6]-EA*alpha*tem
return secf


### 10. Function: Main routine

#### Main routine processing procedure

• Get I / O file name from command line </ b>: Get the input / output data file name as a command line argument. Input data file name: fnameR </ var>, Output data file name: fnameW </ var>

• Data entry </ b>: Get input data by function inpdata_3dfrm </ em>

• Export input data </ b>: Write input data with the function prinp_3dfrm </ em>

• Rigid matrix assembly </ b>: For each element, the stiffness matrix ck </ var> (after coordinate conversion), temperature load vector tfe </ var> (after coordinate conversion), and inertial force vector bfe </ var> Calculate and incorporate into the overall stiffness equation. Here, the arguments of the functions of the stiffness matrix, temperature load vector, and inertial force vector calculation are all scalars. If the element number is specified, the node coordinates and element physical property values that make up the element can be specified, so it seems that sending a scalar is more programmatically cleaner than an array that contains the values of other elements. .. The overall stiffness matrix is added to the variable gk </ var>, and the load vector is added to the variable fp </ var>. The position information for incorporating the element stiffness matrix into the overall stiffness matrix is stored in the array ir [] </ var>.

• Boundary condition processing (processing of known displacement degree) </ b>: The size of the overall stiffness matrix is not changed, and the nodes with known displacement are processed.

Item Description
npoin Total number of nodes less
mpfix Array containing constraint conditions for each node (= 1: constraint, = 0: free)
rdis Displacement amount of the node to be displaced
gk Overall stiffness matrix
disg Solution (displacement) of the overall stiffness equation
• Solving the stiffness equation (displacement calculation) </ b>: After CSR conversion of the overall stiffness matrix gk </ var>, solve simultaneous linear equations with scipy.sparse.linalg.spsolve () </ var>.

• Boundary condition processing (incorporation of known displacement) </ b>: After solving the simultaneous linear equations, remember to re-enter the known displacement in the term of degrees of freedom corresponding to the known displacement.

• Element section force calculation </ b>: Calculation of section force for each element by the function calsecf_3dfrm </ em>

• Export calculation results </ b>: Calculation result output by function prout_3dfrm </ em>

• Information display such as calculation time </ b>: By the following, the difference between the calculation completion time and the start time is output as the processing time.

import time

start=time.time()
......
dtime=time.time()-start


#### Main routine program

def main_3dfrm():
start=time.time()
args = sys.argv
fnameR=args[1] # input data file
fnameW=args[2] # output data file
nod=2              # Number of nodes per element
nfree=6            # Degree of freedom per node
# data input
npoin,nele,nsec,npfix,nlod,ae,node,x,deltaT,mpfix,rdis,fp=inpdata_3dfrm(fnameR,nod,nfree)
# print out of input data
prinp_3dfrm(fnameW,npoin,nele,nsec,npfix,nlod,ae,node,x,deltaT,mpfix,rdis,fp)
# array declaration
ir=np.zeros(nod*nfree,dtype=np.int)   # Work vector for matrix assembly
gk=np.zeros((nfree*npoin,nfree*npoin),dtype=np.float64)   # Global stiffness matrix
# assembly of stiffness matrix & load vectors
for ne in range(0,nele):
i=node[0,ne]-1
j=node[1,ne]-1
m=node[2,ne]-1
x1=x[0,i]; y1=x[1,i]; z1=x[2,i]
x2=x[0,j]; y2=x[1,j]; z2=x[2,j]
ee   =ae[0,m]  # elastic modulus
po   =ae[1,m]  # Poisson's ratio
aa   =ae[2,m]  # section area
aix  =ae[3,m] # tortional constant
aiy  =ae[4,m] # moment of inertia around y-axis
aiz  =ae[5,m] # moment of inertia around z-axis
alpha=ae[7,m] # thermal expansion coefficient
gamma=ae[8,m] # unit weight of material
gkX  =ae[9,m]   # seismic coefficient in X-direction
gkY  =ae[10,m]  # seismic coefficient in Y-direction
gkZ  =ae[11,m]  # seismic coefficient in Z-direction
A=aa  # section area
EA=ee*aa
GJ=ee/2/(1+po)*aix
EIy=ee*aiy
EIz=ee*aiz
tem=0.5*(deltaT[i]+deltaT[j])                # average temperature change
tt   =tm_3dfrm(theta,x1,y1,z1,x2,y2,z2)         # Transformation matrix
ek   =sm_3dfrm(EA,GJ,EIy,EIz,x1,y1,z1,x2,y2,z2) # Stiffness matrix in local coordinate
ck   =np.dot(np.dot(tt.T,ek),tt)             # Stiffness matrix in global coordinate
tfe_l=tfvec_3dfrm(EA,alpha,tem)              # Thermal load vector in local coordinate
tfe  =np.dot(tt.T,tfe_l)                     # Thermal load vector in global coordinate
bfe  =bfvec_3dfrm(A,gamma,gkX,gkY,gkZ,x1,y1,z1,x2,y2,z2) # Body force vector in global coordinate
ir[11]=6*j+5; ir[10]=ir[11]-1; ir[9]=ir[10]-1; ir[8]=ir[9]-1; ir[7]=ir[8]-1; ir[6]=ir[7]-1
ir[5] =6*i+5; ir[4] =ir[5]-1 ; ir[3]=ir[4]-1 ; ir[2]=ir[3]-1; ir[1]=ir[2]-1; ir[0]=ir[1]-1
for i in range(0,nod*nfree):
it=ir[i]
fp[it]=fp[it]+tfe[i]+bfe[i]
for j in range(0,nod*nfree):
jt=ir[j]
gk[it,jt]=gk[it,jt]+ck[i,j]
# treatment of boundary conditions
for i in range(0,npoin):
for j in range(0,nfree):
if mpfix[j,i]==1:
iz=i*nfree+j
fp[iz]=0.0
for i in range(0,npoin):
for j in range(0,nfree):
if mpfix[j,i]==1:
iz=i*nfree+j
fp=fp-rdis[j,i]*gk[:,iz]
gk[:,iz]=0.0
gk[iz,iz]=1.0
# solution of simultaneous linear equations
#disg = np.linalg.solve(gk, fp)
gk = csr_matrix(gk)
disg = spsolve(gk, fp, use_umfpack=True)
# recovery of restricted displacements
for i in range(0,npoin):
for j in range(0,nfree):
if mpfix[j,i]==1:
iz=i*nfree+j
disg[iz]=rdis[j,i]
# calculation of section force
dis=np.zeros(12,dtype=np.float64)
fsec  =np.zeros((nod*nfree,nele),dtype=np.float64)  # Section force vector
for ne in range(0,nele):
i=node[0,ne]-1
j=node[1,ne]-1
m=node[2,ne]-1
x1=x[0,i]; y1=x[1,i]; z1=x[2,i]
x2=x[0,j]; y2=x[1,j]; z2=x[2,j]
ee   =ae[0,m]  # elastic modulus
po   =ae[1,m]  # Poisson's ratio
aa   =ae[2,m]  # section area
aix  =ae[3,m] # tortional constant
aiy  =ae[4,m] # moment of inertia around y-axis
aiz  =ae[5,m] # moment of inertia around z-axis
alpha=ae[7,m] # thermal expansion coefficient
EA=ee*aa
GJ=ee/2/(1+po)*aix
EIy=ee*aiy
EIz=ee*aiz
tem=0.5*(deltaT[i]+deltaT[j]) # average temperature change
dis[0]=disg[6*i]  ; dis[1] =disg[6*i+1]; dis[2]= disg[6*i+2]
dis[3]=disg[6*i+3]; dis[4] =disg[6*i+4]; dis[5]= disg[6*i+5]
dis[6]=disg[6*j]  ; dis[7] =disg[6*j+1]; dis[8]= disg[6*j+2]
dis[9]=disg[6*j+3]; dis[10]=disg[6*j+4]; dis[11]=disg[6*j+5]
fsec[:,ne]=calsecf_3dfrm(EA,GJ,EIy,EIz,theta,alpha,tem,dis,x1,y1,z1,x2,y2,z2)
# print out of result
prout_3dfrm(fnameW,npoin,nele,node,disg,fsec)
# information
dtime=time.time()-start
print('n={0}  time={1:.3f}'.format(nfree*npoin,dtime)+' sec')
fout=open(fnameW,'a')
print('n={0}  time={1:.3f}'.format(nfree*npoin,dtime)+' sec',file=fout)
fout.close()


### 11. Main routine execution

#==============
# Execution
#==============
if __name__ == '__main__': main_3dfrm()


that's all