Python：FEMでの疎行列計算利用による高速化（３次元骨組構造解析での事例）

## １．はじめに
最近，自作の３次元骨組構造解析プログラムで計算を行っているが，自由度が大きくなると計算に時間がかかるのが気になりだしました．
6000自由度（行列サイズは 6000 x 6000）となると，numpy.linalg.solve では計算時間は10秒近くかかります．これを何十ケースもやるとなると，コンピュータの前で待っているのは手持ち無沙汰です．そこで何か早くする方法はないものかと探っていて見つけた方法が，scipyの疎行列計算の活用です．

結果からいうと，下表の通り，疎行列処理の導入により，十分に早くなります．自由度が小さく１秒以下の計算での違いなどは気になりませんが，自由度が大きい場合に１ケース９秒が３秒になる時間短縮は大きい！そこで，この記事では，疎行列計算ライブラリを用いた連立一次方程式解法の簡単な例と，実際のFEMプログラムへの実装を説明します．実装と行っても，プログラムの冒頭でライブラリを読み込んで，少し行列を変形するだけですが．．．

| 自由度 (行列サイズ：n x n) | n = 864 | n = 2166 | n = 6012 | 備考 |
|:------:|:-------:|:--------:|:--------:|:--------:|
| numpy.linalg.solve          | 0.125 sec | 0.677 sec | 8.909 sec | 疎行列処理なし |
| scipy.sparse.linalg.spsolve | 0.123 sec | 0.411 sec | 2.645 sec | 疎行列処理あり |

ちなみにどのくらい疎か，n=6012 の例で見てみると．．．たしかに疎だ．これを縦横全部チェックするのは非効率だし，バンド化しても空白域は大きい．疎行列処理のテクニックを導入するのがよかろうと思うわけです．なお，拘束節点処理後なのでもちろん対称ではありません．もし対称だったら大問題．．．特異行列で解けません．

<img src="https://qiita-image-store.s3.amazonaws.com/0/129300/244759ac-6bc1-7a84-45bd-af9f007a6c98.png" width=70%>

### （関連リンク）

+ [Pythonで立体骨組構造解析](https://qiita.com/damyarou/items/2b83b9b20181c1b8e35d)
+ [Pythonによる立体骨組解析事例](https://qiita.com/damyarou/items/6fddb01f318190967084)
+ [有限要素法とPythonプログラム（三角形一次要素を用いた２次元応力解析プログラム）](https://qiita.com/damyarou/items/8ca3432f9bba20a1f5ea)

## ２．Scipy 疎行列ライブラリの利用
以下のサイトを参考にしました．

+ [http://www.turbare.net/transl/scipy-lecture-notes/advanced/scipy_sparse/solvers.html#sparse-direct-solvers](http://www.turbare.net/transl/scipy-lecture-notes/advanced/scipy_sparse/solvers.html#sparse-direct-solvers)
+ [https://ohke.hateblo.jp/entry/2018/01/07/230000](https://ohke.hateblo.jp/entry/2018/01/07/230000)
+ [http://hamukazu.com/2014/09/26/scipy-sparse-basics/](http://hamukazu.com/2014/09/26/scipy-sparse-basics/)
+ [http://hamukazu.com/2014/12/03/internal-data-structure-scipy-sparse/](http://hamukazu.com/2014/12/03/internal-data-structure-scipy-sparse/)

### 簡単な連立一次方程式を解いてみる（１）
上述サイトを参考に，以下のプログラムを作成して実行．

```py:
import numpy as np
from scipy.sparse.linalg import dsolve
aa=np.array([[ 1,  5,  0,  0,  0],
             [ 0,  2,  8,  0,  0],
             [ 0,  0,  3,  9,  0],
             [ 0,  0,  0,  4, 10],
             [ 0,  0,  0,  0,  5]],dtype=np.float64)
bb=np.array([1, 2, 3, 4, 5], dtype=np.float64)
x = dsolve.spsolve(aa, bb, use_umfpack=True)
print('x=',x)
print('Error=',np.dot(aa,x)-bb)  
```

結果は以下の通り．解けていますが，警告が出ています．「spsplve に放り込むには，CSC か CSR マトリックスのフォーマットにしてね．」ということのようです．私は記法の違いを完全に理解していませんが，とりあえず「CSR」を採用[(こちらが参考になります)](http://hamukazu.com/2014/12/03/internal-data-structure-scipy-sparse/)．

```py:
x= [106.  -21.    5.5  -1.5   1. ]
Error= [0. 0. 0. 0. 0.]

/usr/local/lib/python3.7/site-packages/scipy/sparse/linalg/dsolve/linsolve.py:133: SparseEfficiencyWarning: spsolve requires A be CSC or CSR matrix format
  SparseEfficiencyWarning)
```

### 簡単な連立一次方程式を解いてみる（２）
次に以下を実行．lil_matrix と csr 形式に書き換える tocsr() を使うため，
<em>from scipy.sparse import lil_matrix, csr_matrix</em> をインポート．
lil_matrix は，次元は元の行列と同じですが，非ゼロ項のみを記憶するものらしい．

```py
import numpy as np
from scipy.sparse.linalg import dsolve
from scipy.sparse import lil_matrix, csr_matrix
aa=np.array([[ 1,  5,  0,  0,  0],
             [ 0,  2,  8,  0,  0],
             [ 0,  0,  3,  9,  0],
             [ 0,  0,  0,  4, 10],
             [ 0,  0,  0,  0,  5]],dtype=np.float64)
cc=aa
bb = np.array([1, 2, 3, 4, 5], dtype=np.float64)

aa=lil_matrix(aa).tocsr()
print(aa.shape)
print(aa)
x = dsolve.spsolve(aa, bb, use_umfpack=True)
print('x=',x)
print('Error=',np.dot(cc,x)-bb)  
```

結果は以下の通り．警告なしで正解が出ている模様．

```py
(5, 5)
  (0, 0)	1.0
  (0, 1)	5.0
  (1, 1)	2.0
  (1, 2)	8.0
  (2, 2)	3.0
  (2, 3)	9.0
  (3, 3)	4.0
  (3, 4)	10.0
  (4, 4)	5.0
x= [106.  -21.    5.5  -1.5   1. ]
Error= [-1.42108547e-14  0.00000000e+00 -3.55271368e-15  0.00000000e+00
  0.00000000e+00]
 ```

### 簡単な連立一次方程式を問いてみる（３）
scipy のドキュメントを見ていたら，csr_matrix については，lil_matrix を介さなくてもそのまま変換できそうであること，
spsolve についても dsolve を介さなくても使えそうなことがわかったので，記述量を減らしました．
 
+ scipy.sparse.csr_matrix [https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csr_matrix.html#scipy.sparse.csr_matrix](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csr_matrix.html#scipy.sparse.csr_matrix)  
+ scipy.sparse.linalg.spsolve [https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.spsolve.html#scipy.sparse.linalg.spsolve](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.spsolve.html#scipy.sparse.linalg.spsolve)

```py:
import numpy as np
from scipy.sparse.linalg import spsolve
from scipy.sparse import csr_matrix
aa=np.array([[ 1,  5,  0,  0,  0],
             [ 0,  2,  8,  0,  0],
             [ 0,  0,  3,  9,  0],
             [ 0,  0,  0,  4, 10],
             [ 0,  0,  0,  0,  5]],dtype=np.float64)
cc=aa
bb = np.array([1, 2, 3, 4, 5], dtype=np.float64)

aa=csr_matrix(aa)
print(aa.shape)
print(aa)
x = spsolve(aa, bb, use_umfpack=True)
print('x=',x)
print('Error=',np.dot(cc,x)-bb)  
```
出力はもちろん（２）と同じです．

```py:
(5, 5)
  (0, 0)	1.0
  (0, 1)	5.0
  (1, 1)	2.0
  (1, 2)	8.0
  (2, 2)	3.0
  (2, 3)	9.0
  (3, 3)	4.0
  (3, 4)	10.0
  (4, 4)	5.0
x= [106.  -21.    5.5  -1.5   1. ]
Error= [-1.42108547e-14  0.00000000e+00 -3.55271368e-15  0.00000000e+00
  0.00000000e+00]
```

### 結局どうするのか？　上述（３）をベースに．

（a）以下をインポート

```py
import numpy as np
from scipy.sparse.linalg import spsolve
from scipy.sparse import csr_matrix
```

（b）連立一次方程式  $\\{b\\}=[A]\\{x\\}$  の行列 $aa=[A]$ とベクトル $bb=\\{b\\}$ を numpy 配列で作成．

（c）以下に示すように，行列 $aa$ を <em>csr_matrix</em> で変形し，ベクトル $bb$ とともに <em>dsolve.spsolve</em> に突っ込めば終わり．

```py
aa=csr_matrix(aa)
x = spsolve(aa, bb, use_umfpack=True)
```

### 付録：どのくらい疎なのかを見てみる
CSR 形式の行列を保存するには，

```py:
numpy.save('array',gk) # arrayはファイル名
```

でOK[ (numpy.save) ](https://docs.scipy.org/doc/numpy-1.14.0/reference/generated/numpy.save.html)．array.npy というバイナリファイルが出来上がります．これを読み込んで非ゼロ要素を表示するプログラムは以下の通り．前出の図を描くプログラムです．CSR を toarray でフルの正方行列に戻し，for 文を２重にまわして非ゼロ要素を探しているので，遅いです．もっといい方法がありそうなのだけど，今の所思いつかない．

```py:
# np.save('array',gk): save an array to a binary file in Numpy .npy format. 
# a .npy extension is appended to the file name.
import numpy as np
import matplotlib.pyplot as plt

gk=np.load('array.npy')
gk=gk.any().toarray()

n=gk.shape[0]
m=gk.shape[1]
x=[]
y=[]
for i in range(0,n):
    for j in range(0,m):
        if 0<abs(gk[i,j]):
            x=x+[float(i)]
            y=y+[float(j)]

fsz=12
xmin=0
xmax=float(m-1)
ymin=0
ymax=float(n-1)

fig=plt.figure(figsize=(6,6),facecolor='w')
plt.rcParams['font.size']=fsz
plt.rcParams['font.family']='sans-serif'

plt.xlim([xmin,xmax])
plt.ylim([ymax,ymin])
plt.xlabel('Column')
plt.ylabel('Row')
plt.gca().set_aspect('equal',adjustable='box')
plt.plot(x,y,'s',ms=1)
strt='Size: {0:d} x {1:d}, Non-zero: {2:d}'.format(n,m,len(x))
plt.title(strt,loc='left',fontsize=fsz)

fnameF='fig_sparse.png'
plt.savefig(fnameF, dpi=200, bbox_inches="tight", pad_inches=0.1)
plt.show()
```


## ３．実プログラム（３次元骨組構造解析プログラム）
[こちらの記事](https://qiita.com/damyarou/items/2b83b9b20181c1b8e35d)で紹介した，３次元骨組構造解析プログラムを改造したものを紹介します．昔のプログラムからは，以下の点をいじっています．


+ 疎行列処理導入による高速化
+ 関数名を小文字にした
+ データ入力部を関数にした
+ 計算用関数（データ入力・プリントアウト部を除く）の引数を全てスカラーとした
+ 無意味な配列初期化の部分をカットした
+ 配列初期化時の形状をタプルで指定するようにした（何故かこれまでリストで渡していた）
+ メインルーチンも関数化した．これは下の記事を見てハット思い実行しました．これによりグローバル変数となっているものを関数の引数で呼び込んだりというおかしなことがなくなったと思います．[（pythonで小さなツールを作る時のtips）](https://qiita.com/m_mizutani/items/c48f67f871d1d41ff4b9)

この疎行列を扱う scipy ライブラリを使う方法は，２次元応力解析など，自由度が大きなモデルの連立一次方程式を解くのに有効に使えそうです．

### 疎行列導入に伴い変えるところ

splolve と csr_matrix のインポートを追加．

```py:
from scipy.sparse.linalg import spsolve
from scipy.sparse import csr_matrix
```

剛性マトリックスの境界条件処理後，全体剛性マトリックス gk を CSR フォーマットに書き換えて，spsolve に渡します．numpy の np.linalg.solve はコメントアウトしています．

```py
# solution of simultaneous linear equations
#disg = np.linalg.solve(gk, fp)
gk = csr_matrix(gk)
disg = spsolve(gk, fp, use_umfpack=True)
```

### プログラム全文

```py
# ======================
# 3D Frame Analysis
# ======================
import numpy as np
from scipy.sparse.linalg import spsolve
from scipy.sparse import csr_matrix
import sys
import time


def inpdata_3dfrm(fnameR,nod,nfree):
    f=open(fnameR,'r')
    text=f.readline()
    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=f.readline()
        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=f.readline()
        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=f.readline()
        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=f.readline()
        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
    # load
    if 0<nlod:
        for i in range(0,nlod):
            text=f.readline()
            text=text.strip()
            text=text.split()
            lp=int(text[0])           #loaded node
            fp[6*lp-6]=float(text[1]) #load in x-direction
            fp[6*lp-5]=float(text[2]) #load in y-direction
            fp[6*lp-4]=float(text[3]) #load in z-direction
            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


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:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>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:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s} {7:>15s}'
    .format('sec','E','po','A','J','Iy','Iz','theta'),file=fout)
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>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:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s} {7:>15s} {8:>15s} {9:>15s} {10:>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:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>5s} {5:>5s} {6:>5s} {7:>15s} {8:>15s} {9:>15s} {10:>15s} {11:>15s} {12:>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<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:>5s} {1:>5s} {2:>5s} {3:>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()


def prout_3dfrm(fnameW,npoin,nele,node,disg,fsec):
    fout=open(fnameW,'a')
    # displacement
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>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:>5s} {1:>5s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s} {7:>15s}'
    .format('elem','nodei','N_i','Sy_i','Sz_i','Mx_i','My_i','Mz_i'),file=fout)
    print('{0:>5s} {1:>5s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s} {7:>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()


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


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


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


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


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


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
        theta=np.radians(ae[6,m]) # chord angle
        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
        theta=np.radians(ae[6,m]) # chord angle
        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()


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

以　上