有限要素法(FEM) を python で作成する ～vba→python 翻訳編～


#シリーズ目次

- 有限要素法(FEM) を python で作成する ～vba→python 翻訳編～
- [有限要素法(FEM) を python で作成する ～damyarouさんのいじくる編～](http://qiita.com/sasaco/private/ad06f769353cbe12e313)
- 有限要素法(FEM) を python で作成する ～行列の解法編～（予定）

---

# vba→python 翻訳編 の概要

構想をざっくり説明しますと
有限要素法(FEM) を python で作成しよう

まずは、～vba→python 翻訳編～ としまして[この書籍](https://rpx.a8.net/svt/ejp?a8mat=2TINC2+D4Q7VE+2HOM+BW8O1&rakuten=y&a8ejpredirect=http%3A%2F%2Fhb.afl.rakuten.co.jp%2Fhgc%2F0eac8dc2.9a477d4e.0eac8dc3.0aa56a48%2Fa08100810103_2TINC2_D4Q7VE_2HOM_BW8O1%3Fpc%3Dhttp%253A%252F%252Fbooks.rakuten.co.jp%252Frb%252F1700147%252F%26m%3Dhttp%253A%252F%252Fbooks.rakuten.co.jp%252Frb%252F1700147%252F)のエクセルサンプル `5-1 三次元骨組解析.xls`(以下サンプル) をpython で書き直してみようと思います。

![41N5A3TAQXL._SX332_BO1,204,203,200_[1].jpg](https://qiita-image-store.s3.amazonaws.com/0/142847/5eafb41d-020e-b2fe-475a-992737d64922.jpeg)

[エクセル有限要素法入門 骨組構造解析編 (Excel土木講座)](https://rpx.a8.net/svt/ejp?a8mat=2TINC2+D4Q7VE+2HOM+BW8O1&rakuten=y&a8ejpredirect=http%3A%2F%2Fhb.afl.rakuten.co.jp%2Fhgc%2F0eac8dc2.9a477d4e.0eac8dc3.0aa56a48%2Fa08100810103_2TINC2_D4Q7VE_2HOM_BW8O1%3Fpc%3Dhttp%253A%252F%252Fbooks.rakuten.co.jp%252Frb%252F1700147%252F%26m%3Dhttp%253A%252F%252Fbooks.rakuten.co.jp%252Frb%252F1700147%252F)

翻訳は下記のように [サンプル](https://rpx.a8.net/svt/ejp?a8mat=2TINC2+D4Q7VE+2HOM+BW8O1&rakuten=y&a8ejpredirect=http%3A%2F%2Fhb.afl.rakuten.co.jp%2Fhgc%2F0eac8dc2.9a477d4e.0eac8dc3.0aa56a48%2Fa08100810103_2TINC2_D4Q7VE_2HOM_BW8O1%3Fpc%3Dhttp%253A%252F%252Fbooks.rakuten.co.jp%252Frb%252F1700147%252F%26m%3Dhttp%253A%252F%252Fbooks.rakuten.co.jp%252Frb%252F1700147%252F) → [python](https://github.com/sasaco/FramePython/tree/%231) に翻訳していきます。


```vbnet:5-1三次元骨組解析.xls

Private Sub cmd計算_Click()

    Dim I As Long
    Dim NEQ As Integer

    For I = 1 To 配列上限
        AJCB(I) = 0
    Next

    Call データ入力
    Call SKYマトリックス(NEQ)
    For I = 1 To NEQ
        FORCE(I) = 0
    Next I
    Call 分布荷重振り分け
    Call 外力add
    Call decomp(NEQ)
    Call redbak(NEQ)
    Call 変位計算
    Call 力とモーメントの計算

    Call 結果出力

    MsgBox "解析が終わりました"

End Sub
```

```py:framecalc.py
from Calcrate import Calcrate

if __name__ == "__main__":

    cal = Calcrate()
    NEQ = 0

    cal.データ入力()
    cal.SKYマトリックス(NEQ)
    cal.分布荷重振り分け()
    cal.外力add()
    cal.decomp(NEQ)
    cal.redbak(NEQ)
    cal.変位計算()
    cal.力とモーメントの計算()

    cal.結果出力()

    print("解析が終わりました")
```

## 環境条件
私の環境は下記の通りです。

### OS
- windows 7 Home premium Sp1

### 開発環境
- python 3.5.2
- Anaconda 4.2.0(64-bit)
- pandas 0.18.1
- numpy 1.12.0
- tensorflow 0.12.0re1


## まずは動かしてみよう

はじめにソースコードをダウンロードします。

[サンプル](https://rpx.a8.net/svt/ejp?a8mat=2TINC2+D4Q7VE+2HOM+BW8O1&rakuten=y&a8ejpredirect=http%3A%2F%2Fhb.afl.rakuten.co.jp%2Fhgc%2F0eac8dc2.9a477d4e.0eac8dc3.0aa56a48%2Fa08100810103_2TINC2_D4Q7VE_2HOM_BW8O1%3Fpc%3Dhttp%253A%252F%252Fbooks.rakuten.co.jp%252Frb%252F1700147%252F%26m%3Dhttp%253A%252F%252Fbooks.rakuten.co.jp%252Frb%252F1700147%252F) `5-1 三次元骨組解析.xls` は本の案内にしたがってダウンロードしてください。
このエクセルVBA を python に翻訳していきます。

python に翻訳したソースはここにおいておきます。
[`$ git clone -b "#1" https://github.com/sasaco/FramePython.git`](https://github.com/sasaco/FramePython.git)


環境が整ったら、ソースコードのディレクトリに移動して、framecalc.pyを叩くと解析が始まります。
"解析が終わりました" と出れば、正しく解析が行われました。

```text:
$ cd FramePython

$ python framecalc.py

$ 解析が終わりました
```

# python プログラムの仕様

## データの入出力について

[サンプル](https://rpx.a8.net/svt/ejp?a8mat=2TINC2+D4Q7VE+2HOM+BW8O1&rakuten=y&a8ejpredirect=http%3A%2F%2Fhb.afl.rakuten.co.jp%2Fhgc%2F0eac8dc2.9a477d4e.0eac8dc3.0aa56a48%2Fa08100810103_2TINC2_D4Q7VE_2HOM_BW8O1%3Fpc%3Dhttp%253A%252F%252Fbooks.rakuten.co.jp%252Frb%252F1700147%252F%26m%3Dhttp%253A%252F%252Fbooks.rakuten.co.jp%252Frb%252F1700147%252F)エクセルは シート**[入力1][入力2]**をインプットデータとし,
シート**[計算・出力]**に解析結果を格納しています。
![b.png](https://qiita-image-store.s3.amazonaws.com/0/142847/35115915-f098-a0b3-cdc0-3ce783fc08ca.png)

翻訳した python プログラムは、エクセルのシートと連携してないので、
[サンプル](https://rpx.a8.net/svt/ejp?a8mat=2TINC2+D4Q7VE+2HOM+BW8O1&rakuten=y&a8ejpredirect=http%3A%2F%2Fhb.afl.rakuten.co.jp%2Fhgc%2F0eac8dc2.9a477d4e.0eac8dc3.0aa56a48%2Fa08100810103_2TINC2_D4Q7VE_2HOM_BW8O1%3Fpc%3Dhttp%253A%252F%252Fbooks.rakuten.co.jp%252Frb%252F1700147%252F%26m%3Dhttp%253A%252F%252Fbooks.rakuten.co.jp%252Frb%252F1700147%252F)エクセルの シート**[入力1][入力2][計算・出力]**を csv に変換したファイルを入出力に用います。

- [入力1.csv](https://github.com/sasaco/FramePython/blob/%231/%E5%85%A5%E5%8A%9B1.csv)
- [入力2.csv](https://github.com/sasaco/FramePython/blob/%231/%E5%85%A5%E5%8A%9B2.csv)
- [計算・出力.csv](https://github.com/sasaco/FramePython/blob/%231/%E8%A8%88%E7%AE%97%E3%83%BB%E5%87%BA%E5%8A%9B.csv)

## 配列の添字について

vba の配列は 1～

```vbnet:5-1三次元骨組解析.xls
Range(1,1)
```

python の配列は 0～

```py:Calcrate.py
Range[0][0]
```

なので、翻訳する際に注意が必要だった。


## 描画機能について
[サンプル](https://rpx.a8.net/svt/ejp?a8mat=2TINC2+D4Q7VE+2HOM+BW8O1&rakuten=y&a8ejpredirect=http%3A%2F%2Fhb.afl.rakuten.co.jp%2Fhgc%2F0eac8dc2.9a477d4e.0eac8dc3.0aa56a48%2Fa08100810103_2TINC2_D4Q7VE_2HOM_BW8O1%3Fpc%3Dhttp%253A%252F%252Fbooks.rakuten.co.jp%252Frb%252F1700147%252F%26m%3Dhttp%253A%252F%252Fbooks.rakuten.co.jp%252Frb%252F1700147%252F)エクセルの描画機能は python に移植しない

# 解析結果の比較

ひとまず vba→python の翻訳作業が終わったので、両者の解析結果を比較してみました
ほぼ 0 の小さい値 は誤差があるが、まあまあな再現度ではないでしょうか

## サンプルエクセル の 解析結果（抜粋）

| 変位 |          |           |           |           |           |           |   | REACT |           |          |         |         |           |         |   | 要素 |      |          |         |         |           |         |          |
|------|----------|-----------|-----------|-----------|-----------|-----------|---|-------|-----------|----------|---------|---------|-----------|---------|---|------|------|----------|---------|---------|-----------|---------|----------|
|      | X        | Y         | Z         | angX      | angY      | angZ      |   |       | Rx        | Ry       | Rz      | Mx      | My        | Mz      |   |      | 節点 | Fx       | Fy      | Fz      | Mx        | My      | Mz       |
| 1    | 0.00E+00 | 0.00E+00  | 0.00E+00  | 0.00E+00  | 0.00E+00  | 0.00E+00  |   | 1     | 3.49E+01  | 408.073  | 20.088  | 31.400  | 1.73E-16  | -17.837 |   | 1    | 1    | 408.073  | -34.929 | 20.088  | 1.73E-16  | -31.400 | -17.837  |
| 2    | 2.30E-03 | -2.59E-02 | -1.23E-04 | 1.17E-04  | -6.91E-21 | -2.08E-03 |   | 2     | -5.00E-15 | -76.000  | 0.000   | 48.000  | 1.28E-14  | -75.000 |   |      | 2    | -408.073 | 34.929  | -20.088 | -1.73E-16 | -48.953 | -121.877 |
| 3    | 1.17E-02 | -4.30E-02 | 2.04E-04  | 2.42E-04  | 4.04E-19  | -3.80E-03 |   | 3     | -1.46E-13 | -78.000  | 0.000   | 48.000  | -3.70E-14 | -75.000 |   | 2    | 2    | 269.653  | -44.737 | 22.027  | -1.03E-14 | -34.644 | -46.608  |
| 4    | 2.76E-02 | -5.04E-02 | 1.19E-03  | 4.20E-04  | -5.62E-19 | -5.41E-03 |   | 4     | 9.19E-14  | -56.000  | 0.000   | 42.667  | 3.56E-13  | -66.667 |   |      | 3    | -269.653 | 44.737  | -22.027 | 1.03E-14  | -53.464 | -132.339 |
| 5    | 0.00E+00 | 0.00E+00  | 0.00E+00  | 0.00E+00  | 0.00E+00  | 0.00E+00  |   | 5     | 7.02E+00  | 623.968  | 20.088  | 31.400  | -6.73E-15 | 29.101  |   | 3    | 3    | 116.788  | -45.660 | 18.809  | 2.42E-14  | -24.329 | -50.990  |


## pythonプログラム の 解析結果（抜粋）

| 変位 |             |              |              |              |           |              |   | REACT |              |             |              |              |           |              |   | 要素 |      |              |              |              |           |              |              |
|------|-------------|--------------|--------------|--------------|-----------|--------------|---|-------|--------------|-------------|--------------|--------------|-----------|--------------|---|------|------|--------------|--------------|--------------|-----------|--------------|--------------|
|      | X           | Y            | Z            | angX         | angY      | angZ         |   |       | Rx           | Ry          | Rz           | Mx           | My        | Mz           |   |      | 節点 | Fx           | Fy           | Fz           | Mx        | My           | Mz           |
| 1    | 0           | 0            | 0            | 0            | 0         | 0            |   | 1     | 34.9286048   | 408.0727049 | 20.08821182  | 31.40032555  | -9.49E-15 | -17.83716034 |   | 1    | 1    | 408.0727049  | -34.9286048  | 20.08821182  | -9.49E-15 | -31.40032555 | -17.83716034 |
| 2    | 0.002298747 | -0.025909378 | -0.000123095 | 0.000117015  | 3.79E-19  | -0.002080804 |   | 2     | -1.91E-14    | -76         | -3.11E-15    | 48           | 1.46E-14  | -75          |   |      | 2    | -408.0727049 | 34.9286048   | -20.08821182 | 9.49E-15  | -48.95252174 | -121.8772589 |
| 3    | 0.011665227 | -0.043030177 | 0.000204306  | 0.000242481  | 8.13E-19  | -0.003795417 |   | 3     | -6.95E-14    | -78         | -6.66E-15    | 48           | -1.17E-13 | -75          |   | 2    | 2    | 269.6525828  | -44.73675837 | 22.02695187  | -1.08E-14 | -34.64397073 | -46.60823232 |
| 4    | 0.027638145 | -0.050445274 | 0.001194231  | 0.000419675  | 6.78E-19  | -0.005408655 |   | 4     | -2.39E-13    | -56         | 3.55E-14     | 42.66666667  | -2.39E-13 | -66.66666667 |   |      | 3    | -269.6525828 | 44.73675837  | -22.02695187 | 1.08E-14  | -53.46383676 | -132.3388011 |
| 5    | 0           | 0            | 0            | 0            | 0         | 0            |   | 5     | 7.022385521  | 623.9683925 | 20.08821182  | 31.40032555  | -4.93E-15 | 29.10143338  |   | 3    | 3    | 116.7877789  | -45.66043605 | 18.80913902  | 3.39E-15  | -24.32874619 | -50.9899696  |


# 次回

とりあえず似たような結果が出るところまでは、出来ました。

次回は、
- 解析部分のクラスを分ける
- 行列の解法を python らしく

など考えています。

この記事は 現在進行中の企画ですので、
まだ 完成していません。

現在の プログラムの状態は master ブランチにあります。
https://github.com/sasaco/FramePython/tree/master
上記のうち何かの修正を開始しているかもしれません

興味のある方 はご連絡ください。