Dynamoで地震応答解析

# 目次
[1. はじめに](#はじめに)
[2. ライブラリのインストール](#ライブラリのインストール)
[3. 地震応答解析](#地震応答解析)
[4. グラフの出力](#グラフの出力)
[5. おわりに](#おわりに)



# はじめに
構造系ではGrasshopperが人気ですが、Dynamoでも簡単な応答解析ができないか試してみます。



# ライブラリのインストール
Pythonで使う外部ライブラリをインストールします。

__1. はじめに Dynamo を起動し、下記のコードで CPython のバージョンを確認します。__

```Python3:PythonScript
import sys
OUT = sys.version
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2059946/eb4655eb-b6ca-83f3-f136-ec470d453a20.png)

__2. Anaconda をダウンロードし、Anaconda Navigator から「CMD.exe Prompt」を実行します。__

https://www.anaconda.com/products/individual

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2059946/e065b2a0-6bae-36e5-74b4-36446585bb3c.png)

__3. 下記のコードを入力して Dynamo 用の環境をつくります。__
__3.8.3__ は最初に確認した CPython のバージョンの番号です。

```Python3:Prompt
conda create --name Dynamo383 python=3.8.3
```
__4. Proceed ([y]/n)?と聞かれるので Y を入力します。__
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2059946/b598d5ba-fbd4-6516-af3a-90e36888493e.png)

__5. いま作成した環境に切り替えます。__

```Python3:Prompt
conda activate Dynamo383
```

__6. 必要なライブラリをまとめてインストールします。__

```Python3:Prompt
conda install numpy pandas scipy Pillow matplotlib
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2059946/d2be037b-adb0-afa3-1e02-fa02e3f2c29b.png)



# 地震応答解析
せん断変形のみを考慮した多質点系モデルの振動解析を実装します。
波形データは[東京都建設局](https://www.kensetsu.metro.tokyo.lg.jp/jigyo/tech/start/03-jyouhou/kyousin/hakeikiroku/201103jishinhdata.html)で公開されているCSVデータを利用します。

### 条件値と波形データの入力
ライブラリをインポートします。

```Python3:PythonScript
#ライブラリをロード
sys.path.append(r'C:\Users\shota\.conda\envs\Dynamo383\Lib\site-packages')
import sys, math, io
import numpy as np
import pandas as pd
from scipy import interpolate
import matplotlib.pyplot as plt
from PIL import Image as Img
import System.Drawing
from System.Drawing import *
from System.Drawing.Imaging import *
from System.IO import MemoryStream
```

コード中に __IN[0]__ とかあると分かりにくいので入力値は __Inputクラス__ にまとめておきます。

```Python3:PythonScript
# インプット
class Input:
    def __init__(self):
        self.g = IN[0] # 重力加速度
        self.dt = IN[1] # 積分時間刻み
        self.beta = IN[2] # Newmarkβ法のβ
        self.damp_factor = IN[3] # 減衰定数
        self.natural_frequency = IN[4] # 振動数
        self.wave_path = IN[7] # CSVファイルのパス
        self.arr_m = IN[5] # 重量
        self.arr_k = IN[6] # 剛性
```

CSV形式の波形データを読み込みます。

```Python3:PythonScript
# 波形データ
def WaveData(input):
    with open(input.wave_path, 'r') as f:
        df = pd.read_csv(f)
    return list(df['acc'].values.tolist())
```

### マトリクス生成
質量マトリクス、減衰マトリクス、剛性マトリクスをつくります。

```Python3:PythonScript
# 質量マトリクス
def MassMatrix(arr_m, input):    
    return np.diag(list(map(lambda x: x / input.g, arr_m)))
    
# 剛性マトリクス
def StiffMatrix(arr_k):
    size = len(arr_k)
    full = np.zeros((size, size))
    for i in range(size):
        k = arr_k[i]
        part = [[k , -k], [-k,   k]]
        for i_row in range(len(part)):
            for i_col in range(len(part)):
                target_row = i-1 + i_row
                target_col = i-1 + i_col
                if 0 <= target_row and 0 <= target_col:
                    full[target_row][target_col] += part[i_row][i_col]
    return full
    
# 減衰マトリクス
def DampMatrix(k, input):
    omega = 2 * math.pi * input.natural_frequency
    return 2 * input.damp_factor * k / omega
```

### 質点系モデル作成
前項で生成したマトリクスを用いて質点系モデルをつくります。

```Python3:PythonScript
# 質点系モデル
class Model:
    def __init__(self,input):
        arr_m = input.arr_m
        arr_k = input.arr_k
        self.size = len(arr_m)
        arr_k.reverse()
        arr_m.reverse()
        self.m = MassMatrix(arr_m, input)
        self.k = StiffMatrix(arr_k)
        self.c = DampMatrix(self.k, input)
```

### 解析
Newmarkβ法を用いて数値解析を行い、入力波形のステップごとに加速度・速度・変位を求めます。

```Python3:PythonScript
def dynamic_analysis(model, input):
# 初期化
    m = model.m
    k = model.k
    c = model.c
    acc0 = WaveData(input)    
    unit_vector = np.ones(model.size) 
    pre_acc0 = 0.0
    time = 0.0
    dis = np.zeros(model.size)
    vel = np.zeros(model.size)
    acc = np.zeros(model.size)
    ddis = np.zeros(model.size)
    dvel = np.zeros(model.size)
    dacc = np.zeros(model.size)
    dis_his = {}
    vel_his = {}
    acc_his = {}
    for i in range(0, model.size):
        dis_his[i] = []
        vel_his[i] = []
        acc_his[i] = []
    time_his = []
    # Newmarkβ法で解析
    for i in range(0, len(acc0)):
        kbar = k + (1.0/(2.0*input.beta*input.dt)) * c + (1.0/(input.beta*input.dt**2.0)) * m
        dp1 = -1.0 * m.dot(unit_vector) * (acc0[i] - pre_acc0)
        dp2 = m.dot((1.0/(input.beta*input.dt))*vel + (1.0/(2.0*input.beta))*acc)
        dp3 = c.dot((1.0/(2.0*input.beta))*vel + (1.0/(4.0*input.beta)-1.0)*acc*input.dt)
        dp = dp1 + dp2 + dp3
        ddis = np.linalg.inv(kbar).dot(dp)
        dvel = (1.0/(2.0*input.beta*input.dt))*ddis - (1.0/(2.0*input.beta))*vel - ((1.0/(4.0*input.beta)-1.0))*acc*input.dt
        dacc = (1.0/(input.beta*input.dt**2.0))*ddis - (1.0/(input.beta*input.dt))*vel - (1.0/(2.0*input.beta))*acc
        dis += ddis
        vel += dvel
        acc += dacc
        acc_abs = acc + [acc0[i] for n in range(1,model.size)]
        [dis_his[i].append(x) for i, x in enumerate(dis)]
        [vel_his[i].append(x) for i, x in enumerate(vel)]
        [acc_his[i].append(x) for i, x in enumerate(acc_abs)]
        time_his.append(time)
        time += input.dt
        pre_acc0 = acc0[i]
    return time_his, dis_his, vel_his, acc_his
```

### 結果出力
前項で作成した__dynamic_analysis__を実行して結果を求めます。

```Python3:PythonScript
# メイン
input = Input()
model = Model(input)
time_his, dis_his, vel_his, acc_his = dynamic_analysis(model, input)
```



# グラフの出力
__matplotlib__ でグラフを作成します。
グラフを __Watch Image__ ノードで表示するために、__System.Drawing.Bitmap__ に変換します。

```Python3:PythonScript
# 画像変換
def convertToBitmap(fig):
    rgba_buf = fig.canvas.buffer_rgba()
    (w,h) = fig.canvas.get_width_height()
    npImgArray = np.frombuffer(rgba_buf, dtype=np.uint8).reshape((h,w,4))   
    bitmap_ = None
    # alphaチャンネルを削除
    if npImgArray.ndim == 3 and npImgArray.shape[-1] == 4:
        npImgArray = npImgArray[:, :, :-1]
    img = Img.fromarray(npImgArray, "RGB")
    # 画像をバイナリデータに変換し、メモリ上に保存
    byteIO = io.BytesIO()
    img.save(byteIO, format='BMP')
    byteArr = byteIO.getvalue()
    # convert to Net ByteArray
    netBytes = System.Array[System.Byte](byteArr)
    with MemoryStream(netBytes) as ms:
        bitmap_ = Bitmap(ms)
    return bitmap_
    
# グラフ
fig = plt.figure(figsize=(10,10))

ax1 = fig.add_subplot(3,1,1)
ax1.set_ylabel("acc") # 縦軸のタイトル
ax1.plot(time_his, acc_his[1], label="1F", color="darkturquoise") # 1階のグラフ
ax1.plot(time_his, acc_his[0], label="2F", color='limegreen') # 2階のグラフ
ax1.legend(loc='lower right') # 凡例の位置
ax1.grid(True) # グリッド

ax2 = fig.add_subplot(3,1,2)
ax2.set_ylabel("vel") # 縦軸のタイトル
ax2.plot(time_his, vel_his[1], label="1F", color="darkturquoise") # 1階のグラフ
ax2.plot(time_his, vel_his[0], label="2F", color='limegreen') # 2階のグラフ
ax2.legend(loc='lower right') # 凡例の位置
ax2.grid(True) # グリッド

ax3 = fig.add_subplot(3,1,3)
ax3.set_xlabel("time") # 横軸のタイトル
ax3.set_ylabel("dis") # 縦軸のタイトル
ax3.plot(time_his, dis_his[1], label="1F", color="darkturquoise") # 1階のグラフ
ax3.plot(time_his, dis_his[0], label="2F", color='limegreen') # 2階のグラフ
ax3.legend(loc='lower right') # 凡例の位置
ax3.grid(True) # グリッド

fig.canvas.draw()
OUT = convertToBitmap(fig)
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2059946/5c95e279-c4a1-b015-67f8-9ab10d80d21b.png)

__DateTime.Now__ ノードを使ってアニメーションにするのも面白いと思います。
![animation (1).gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2059946/de98f02f-0ae8-1e6c-923b-0b292b8805f9.gif)



# おわりに
スクリプト単体では使い道がないですが、データ操作やアニメーション出力の箇所はいろいろと応用できそうな気がします。