平面トラスの有限要素法のプログラムをPythonで組んでみました。以前に自分がつくったコードを引っ張り出してきて、Pythonに移植しました。このくらいのボリュームのコードにじっくり取り組むと面白いです。
トラスのプログラムを実際に使うことは少ないかもしれませんが、最初のFEMコードはやはり平面トラスからです。平面トラスが組めれば、有限要素法の計算の流れがわかるので、立体フレームなどへの改良は簡単です。
０． FEMなど数値解析シリーズ
この記事は「FEMなど数値解析シリーズ」の記事です。一連の記事は、以下のリンク集を参照してください。
１．サンプルコードについて
以下、サンプルコードの書式について概要を説明したいと思います。
(1) 可読性
コードの可読性を高めるために、def関数へのデータの流れをできるだけ分かりやすくする必要があります。長いコードはコーディングした直後は読めても、１ヶ月もたつと「なんだっけ？」となります。トラスからフレームのコードも組みたいので、ここでの手抜きはあとにひびきます。
さて、有限要素法の古い本はFORTRANで書かれたものが多いので、図書館でFORTRANのコードと格闘することがあります。FORTRANは文法がシンプルなので気長に読めばなんとかなりますが、まるで古文書を解読しているような気分になれます。
FORTRANでは「COMMON」「 GLOBAL」を使う方法もありますが、基本的に引数で渡さなければ関数の外の変数を参照できません。また、引数で渡した変数は関数の外と中で同じ変数を示すようになるため、関数内での変更が関数の外に影響を与えます。
プログラム言語によってデータの流れが微妙に違うので、移植する場合にはどのようにすれば可読性が高くなるのか試行錯誤することになります。
Pythonでは引数なしで外側の変数を内側から参照できる上に、global変数を多用すれば関数の引数をなしにすることも容易です。でも、データの流れはわかりにくくなると思います。
そこで、私は以下のような方針でコードを記述することにしました。
- global変数は不用意には使わない。
- 文法上不要であっても、関数に渡すデータは引数に記述する。例えばサンプルコードのdef k_elem(iii, indv, xx, eyg, are)では、iii, indv, xx, eyg, areはglobal変数なので引数は文法上不要ですが、関数の入力を明確にするために記述しました。
- 関数内で変更した変数はできるだけreturanで返す。Pythonでは2つ以上の変数でもreturanで返すことができるので、かなり柔軟に対応できます。
- やむを得ずglobal変数を用いる場合は、文法上不要であっても関数内でglobal宣言を行う。
- 引数で入力した変数の変更が関数の外に影響を与える場合は、文法上不要であってもその変数をreturnで返す。例えば関数set_stiff()のgkは、引数にgkを記述するだけでも文法的には問題ないが、returnでgkを返して計算の出力の範囲を明確にする。
なお、サンプルコードで上記の方針から外れるのは以下の通りです。
- 入力データを扱うdef datum()、計算の準備データを扱うbc_index()は、あまりにも関数から返すデータが多いため、global変数としてデータを関数の外に渡す。
- データの出力を行うoutput_datum()、output_result()も、あまりにも関数に渡すデータが多いため、引数では渡さずglobal変数としてデータを関数に渡す。 また、ただの出力であるためglobal宣言も行わない。
(2) その他の書式についての備考
- あとで読み返すときに時間がかかるのが嫌いなのでコメント文は多めです。多分、多すぎると思います。本来は解説資料を別に作り、説明はシンプルにする方がプログラムの全体が見渡せてよいと思います。
- プログラムを使うにはファイルで入出力を行った方が楽です。でも、このトラスプログラムの用途はFEMの勉強ですので、小回りがきくようにデータを変数に直接代入しました。
- 少数（float）の出力書式で桁数の指定は行っていません。これはコードを試行錯誤する上で小数点以下の誤差がわかるようにしたかったからです。
- 最も解析に時間がかかるのは連立一次方程式の部分です。このサンプルコードではnumpy.linalg.solve()を使っています。NumPyは計算速度に定評のあるライブラリであり、たった１行で解が求まるのは本当に楽です。なお、規模の大きい構造物など解析時間が気になる場合には、反復法や、疎行列であることを利用した計算方法などを考えるとよいです。
- 入力データ例のdef datum()、def datum_exp01()において、ヤング係数eygは計算内容の検証で値がシンプルになるように1.0を入力しています。
２．サンプルコード
サンプルコードの紹介です。なお、コードの解説は、別の記事で追加していきたいと思います。
# -*- coding: utf-8 -*-
'''平面トラス解析プログラム Ver.1.03 (2020.3.11) by Snow Tree in June
平面トラスの節点変位、部材応力、反力の算定
'''
import numpy as np
def datum():
'''データ入力
[global変数として入力データを受け渡し]
ndim(=2):ディメンジョン。自由度はx方向、y方向の2種類
uint_max(=65535):'uint16' 0～65535
nod:総節点数
nel:総部材数
nfix:境界条件を与える節点の数
nload:外力が作用する節点の数
eyg :ヤング係数
xx[i, j]:i節点の座標(j=0:x方向, j=1:y方向)
are[i] :i部材の断面積
indv[i, j]:i部材を構成する節点番号(j=0, 1)
npfix[i]:境界条件を与えられている節点番号
ndfix[i, j]:fixの方向:0, freeの方向:1(j=0:x方向, J=1:y方向)
npload[i]:外力が作用する節点番号
aload[i, j]:外力の大きさ(j=0:x方向, J=1:y方向)
'''
global ndim, uint_max
global nod, nel, nfix, nload
global eyg
global xx, are, indv, npfix, ndfix, npload, aload
ndim = 2
uint_max = 65535
nod = 3
nel = 3
nfix = 2
nload = 1
eyg = 1.
#'uint16' 0~65535
#numpy配列の宣言
xx = np.zeros((nod, ndim), dtype = 'float64')
are = np.zeros(nel, dtype = 'float64')
indv = np.zeros((nel, 2), dtype = 'uint16')
npfix = np.zeros(nfix, dtype = 'uint16')
ndfix = np.zeros((nfix, ndim), dtype = 'bool')
npload = np.zeros(nload, dtype = 'uint16')
aload = np.zeros((nload, ndim), dtype = 'float64')
xx[:] = [[-100, 0],
[0, 100],
[100, 0]]
are[:] = 100.
indv[:] = [[0,1],
[1,2],
[0,2]]
npfix[:] = [2,0]
ndfix[:] = [[1, 0],
[0, 0]]
npload[0] = 1
aload[0, :] = [0., -3.]
def datum_exp01():
'''データ入力例
5節点、7部材トラスの例
'''
global ndim, uint_max
global nod, nel, nfix, nload
global eyg
global xx, are, indv, npfix, ndfix, npload, aload
ndim = 2
uint_max = 65535
nod = 5
nel = 7
nfix = 2
nload = 2
eyg = 1.
#'uint16' 0~65535
#numpy配列の宣言
xx = np.zeros((nod, ndim), dtype = 'float64')
are = np.zeros(nel, dtype = 'float64')
indv = np.zeros((nel, ndim), dtype = 'uint16')
npfix = np.zeros(nfix, dtype = 'uint16')
ndfix = np.zeros((nfix, ndim), dtype = 'bool')
npload = np.zeros(nload, dtype = 'uint16')
aload = np.zeros((nload, ndim), dtype = 'float64')
xx[:] = [[0., 0.],
[100., 100.],
[200., 0.],
[300., 100.],
[400., 0.]]
are[:] = 100.
indv[:] = [[0, 1],
[1, 2],
[2, 3],
[3, 4],
[0, 2],
[2, 4],
[1, 3]]
npfix = [0, 4]
ndfix[:] = [[0, 0],
[1, 0]]
npload[:] = [1, 3]
aload[:] = [[0., -3.],
[0., -3.]]
def output_datum():
'''入力データの出力
'''
print('========================================')
print('１．入力データ')
print('(1) 制御データ')
print(' ndim, uint_max：[{}, {}]'.format(ndim, uint_max))
print('(2) 共通データ')
print(' E：[{}]'.format(eyg))
print('(3) 節点座標')
print(' 節点数：{}'.format(nod))
for i in range(nod):
print('{:4d}: [x:{:12.5f}, y:{:12.5f}]'\
.format(i, xx[i, 0], xx[i, 1]))
print('(4) 部材')
print(' 部材数：{}'.format(nel))
for i in range(nel):
print('{:4d}: (node:{:4d} -- {:4d}), [A:{:12.5f}]'\
.format(i, indv[i, 0], indv[i, 1], are[i]))
print('(5) 境界条件')
print(' 境界条件を与える節点の数：{}'.format(nfix))
print(' [境界条件：fix:0, free:1]')
for i in range(nfix):
print('{:4d}: (node:{:4d}), [{:2d}, {:2d}]'\
.format(i, npfix[i], ndfix[i, 0], ndfix[i, 1]))
print('(6) 外力')
print(' 外力が作用する節点の数：{}'.format(nload))
for i in range(nload):
print('{:4d}: (node:{:4d}), [x:{:12.5f}, y:{:12.5f}]'\
.format(i, npload[i], aload[i, 0], aload[i, 1]))
def bc_index(ndim, uint_max, nod, nfix, npfix, ndfix):
'''計算準備データの生成
[global変数として結果を受け渡し]
nmat:節点数(nod) * ディメンジョン(ndim)
n_free:境界条件がfreeの部分の配列の大きさ
n_fix:境界条件がfixの部分の配列の大きさ
ind_free[nmat → n_free]:free部分の通し番号と節点番号の対応表
ind_fix[n_fix]:fix部分の通し番号と節点番号の対応表
'''
global nmat, n_free, n_fix, ind_free, ind_fix
nmat = nod * ndim
if nmat >= uint_max:
print('節点数が制限値を超えている')
n_fix = 0
for i in range(nfix):
for j in range(ndim):
if ndfix[i, j] == 0:
n_fix += 1
n_free = nmat - n_fix
ind_free = np.arange(nmat, dtype = 'uint16')
ind_fix = np.zeros(n_fix, dtype = 'uint16')
#fixの部分をuint_maxに入れ替え
for i in range(nfix):
ii = npfix[i]
for j in range(ndim):
if ndfix[i, j] == 0:
ind_free[ndim * ii + j] = uint_max
#fixの部分を除き、前につめる。[0～mm]の範囲が必要なデータ
s0 = 0
s1 = 0
for i in range(nmat):
if ind_free[i] == uint_max:
ind_fix[s0] = i
s0 += 1
else:
ind_free[s1] = ind_free[i]
s1 += 1
#本来は以下の1行は不要。ind_freeの内容を読みやすくするために削除
ind_free = np.delete(ind_free, slice(n_free, nmat), 0)
def k_elem(iii, indv, xx, eyg, are):
'''iii番目の要素の全体座標系要素剛性マトリクスの生成
[return]
ek_global[4, 4]:全体座標系要素剛性マトリクス, dtype = 'float64'
'''
ek_local = np.zeros((2, 2), dtype = 'float64')
xl, rotation = transform(iii, indv, xx)
#部材座標系要素剛性マトリクス
ek_local[0, 0] = eyg * are[iii] / xl
ek_local[0, 1] = - ek_local[0, 0]
ek_local[1, 0] = ek_local[0, 1]
ek_local[1, 1] = ek_local[0, 0]
#全体座標系要素剛性マトリクス
ek_global = np.dot(rotation.T, np.dot(ek_local, rotation))
return ek_global
def transform(iii, indv, xx):
'''indv[iii]の部材長、[[c,s,0,0],[0,0,c,s]]の算定
[return]
xl: 部材長
rotation:[[cos(θ), sin(θ), 0, 0],
[0, 0, cos(θ), sin(θ)]]
（全体座標系における部材の傾きをθとする。反時計回りを正。）
'''
rotation = np.zeros((2,4), dtype ='float64')
k0 = indv[iii, 0]
k1 = indv[iii, 1]
x21 = xx[k1, 0] - xx[k0, 0]
y21 = xx[k1, 1] - xx[k0, 1]
xl = np.sqrt(x21 * x21 + y21 * y21)
c = x21 / xl
s = y21 / xl
#rotation[:] = [[c, s, 0., 0.], [0., 0., c, s]]
rotation[0, 0] = c
rotation[1, 2] = c
rotation[0, 1] = s
rotation[1, 3] = s
return xl, rotation
def set_stiff(iii, ek_global, indv, gk):
'''要素剛性マトリクスから全体剛性マトリクスの生成（fix,freeを含む）
[return]
gk[nmat, nmat]:全体剛性マトリクス（fix,freeを含む）
'''
m = 2 #0端、1端
for i in range(m):
i_t = ndim * indv[iii, i]
i_e = ndim * i
for j in range(m):
j_t = ndim * indv[iii, j]
j_e = ndim * j
for ij in range(ndim):
ij_t = i_t + ij
ij_e = i_e + ij
for ji in range(ndim):
ji_t = j_t + ji
ji_e = j_e + ji
gk[ij_t, ji_t] += ek_global[ij_e, ji_e]
return gk
def set_load(ndim, nmat, nload, npload, aload):
'''荷重ベクトルの生成（fix,freeを含む）
[return]
f_load[nmat]:荷重ベクトル
'''
f_load = np.zeros(nmat, dtype = 'float64')
if nload != 0:
for i in range(nload):
k = ndim * npload[i]
for j in range(ndim):
f_load[k + j] += aload[i, j]
return f_load
def extract_fix(nmat, n_free, n_fix, ind_free, ind_fix, gk):
'''全体剛性マトリクスgkのfix部分をk_yfixに代入（反力算定用）
[return]
k_yfix[n_fix, n_free]:全体剛性マトリクスgk[nmat, nmat]のfix部分
'''
k_yfix = np.zeros((n_fix, n_free), dtype = 'float64')
for i in range(n_fix):
ii = ind_fix[i]
for j in range(n_free):
jj = ind_free[j]
k_yfix[i, j] = gk[ii, jj]
return k_yfix
def del_stiff_fix(nmat, n_free, ind_free ,f_load, gk):
'''f_load[], gk[,]のfix部分を削除
[return]
f_load[nmat → n_free]:荷重ベクトル(fix部削除)
gk[nmat → n_free, nmat → n_free]:
全体座標系剛性マトリクス(fix部削除)
'''
for i in range(n_free):
ii = ind_free[i]
f_load[i] = f_load[ii]
for j in range(n_free):
jj = ind_free[j]
gk[i, j] = gk[ii, jj]
#行列計算をするために削除。※numpy配列の形状変更は計算時間を要する
f_load = np.delete(f_load, slice(n_free, nmat), 0)
gk = np.delete(gk, slice(n_free, nmat), 0)
gk = np.delete(gk, slice(n_free, nmat), 1)
return f_load, gk
def calc_disp(ndim, nod, n_free, ind_free, y_free):
'''節点変位の計算
[return]
disp_c[i, j]:i節点のj方向の変位
(j=0:x方向変位, J=1:y方向変位, j=2:変位の大きさ)
'''
disp_c = np.zeros((nod, ndim + 1), dtype = 'float64')
for i in range(n_free):
ii = ind_free[i] // ndim
jj = ind_free[i] % ndim
disp_c[ii, jj] = y_free[i]
for i in range(nod):
dis_x = disp_c[i, 0]
dis_y = disp_c[i, 1]
disp_c[i, 2] = np.sqrt(dis_x ** 2 + dis_y ** 2)
return disp_c
def calc_stress(nel, ndim, indv, xx, disp_c):
'''部材軸力の計算
[return]
stress_c[i]:i部材の軸力（引張側正）
'''
em_tmp = np.zeros(4, dtype = 'float64')
stress_c = np.zeros(nel, dtype = 'float64')
for iii in range(nel):
xl, rotation = transform(iii, indv, xx)
for i in range(2):
k1 = indv[iii, i]
for j in range(ndim):
k2 = ndim * i + j
em_tmp[k2] = disp_c[k1, j]
xs = np.dot(rotation, em_tmp)
stress_c[iii] = eyg * are[iii] * (xs[1] - xs[0]) / xl
return stress_c
def calc_react_f(n_fix, ind_fix, k_yfix, y_free):
'''反力の計算
[return]
rf_c[i, j]:(0 ≦ i ＜ n_fix)
(j=0:反力の節点番号, j=1:反力の向き(x:0, y:1), j=2:反力)
'''
rf_c = np.dot(k_yfix, y_free)
rf_dere = [[], []]
for i in range(n_fix):
rf_dere[0].append(ind_fix[i] // ndim)
if (ind_fix[i] % ndim) == 0:
rf_dere[1].append('x')
else:
rf_dere[1].append('y')
return rf_c, rf_dere
def output_result():
'''計算結果の出力
誤差確認のためfloatの出力は書式指定せず
'''
print('========================================')
print('２．計算結果')
print('(1) 節点変位')
print(' 節点: 変位[x:], [y:], [xy:]')
for i in range(nod):
print('{:4d}: [x:{}], [y:{}], [xy:{}]'\
.format(i, disp_c[i, 0], disp_c[i, 1], disp_c[i, 2]))
print('(2) 部材応力')
print(' 部材: (節点0--節点1), 部材応力')
for i in range(nel):
print('{:4d}: ({:4d} -- {:4d}), {}'\
.format(i, indv[i, 0], indv[i, 1], stress_c[i]))
print('(3) 反力')
print(' No.:(節点-方向), 反力')
for i in range(n_fix):
print('{:4d}:({:4d} - {}), {}'\
.format(i, rf_dere[0][i], rf_dere[1][i], rf_c[i]))
if __name__ == '__main__':
'''メインプログラム
[入力、計算準備データ]
datum(), bc_index() にて
[主な計算用データ]
gk[nmat→nfree]:全体剛性マトリクス
f_load[nmat→nfree]:荷重ベクトル
ek_local[2, 2]:要素剛性マトリクス（部材座標系）
ek_global[4, 4]：要素剛性マトリクス（全体座標系）
k_yfix[n_fix, n_free]:変位がfixされた部分の剛性マトリクスの抜き出し
y_free[n_free]：節点変位のベクトル（freeの部分）
[計算結果]
disp_c[nod]:節点変位のベクトル（free,fix）
stress_c[nel]：部材応力
rf_c[n_fix]：反力
rf_dere[n_fix][i]：反力の節点番号と方向(i=0:節点番号、i=1:方向) ※リスト
'''
#データ入力
datum()
#datum_exp01()
#入力データの出力
output_datum()
#計算準備データ
#global nmat, n_free, n_fix, ind_free, ind_fix
bc_index(ndim, uint_max, nod, nfix, npfix, ndfix)
gk = np.zeros((nmat, nmat), dtype = 'float64')
for iii in range(nel):
#全体座標系要素剛性マトリクスの生成
ek_global = k_elem(iii, indv, xx, eyg, are)
#要素剛性マトリクスより全体剛性マトリクスを生成
gk = set_stiff(iii, ek_global, indv, gk)
#荷重ベクトルの生成
f_load = set_load(ndim, nmat, nload, npload, aload)
#fixされているところの剛性マトリクスのみを抽出（反力算定用）
k_yfix = extract_fix(nmat, n_free, n_fix, ind_free, ind_fix, gk)
#f_load[], gk[,]のfixの部分を削除
f_load, gk = del_stiff_fix(nmat, n_free, ind_free, f_load, gk)
#連立一次方程式を解き、拘束条件のない部分の変位を求める
y_free = np.linalg.solve(gk, f_load)
#節点変位の計算
disp_c = calc_disp(ndim, nod, n_free, ind_free, y_free)
#部材応力の計算
stress_c = calc_stress(nel, ndim, indv, xx, disp_c)
#反力の計算
rf_c, rf_dere =calc_react_f(n_fix, ind_fix, k_yfix, y_free)
#計算結果の出力
output_result()３．解析モデルと解析結果
平面トラス解析プログラムのサンプルコードでは、コードの中でdef datum()、def datum_exp01()の2種類の解析モデルを紹介しています。
(1) 解析モデル１
def datum()に示す解析モデルと解析結果です。
以下、出力結果です。節点変位の「xy:」とはx方向変位とy方向変位から求めた合成変位（二乗和平方根）です。
#出力02
========================================
１．入力データ
(1) 制御データ
ndim, uint_max：[2, 65535]
(2) 共通データ
E：[1.0]
(3) 節点座標
節点数：3
0: [x: -100.00000, y: 0.00000]
1: [x: 0.00000, y: 100.00000]
2: [x: 100.00000, y: 0.00000]
(4) 部材
部材数：3
0: (node: 0 -- 1), [A: 100.00000]
1: (node: 1 -- 2), [A: 100.00000]
2: (node: 0 -- 2), [A: 100.00000]
(5) 境界条件
境界条件を与える節点の数：2
[境界条件：fix:0, free:1]
0: (node: 2), [ 1, 0]
1: (node: 0), [ 0, 0]
(6) 外力
外力が作用する節点の数：1
0: (node: 1), [x: 0.00000, y: -3.00000]
========================================
２．計算結果
(1) 節点変位
節点: 変位[x:], [y:], [xy:]
0: [x:0.0], [y:0.0], [xy:0.0]
1: [x:1.5], [y:-5.7426406871192865], [xy:5.935311454452738]
2: [x:3.0], [y:0.0], [xy:3.0]
(2) 部材応力
部材: (節点0--節点1), 部材応力
0: ( 0 -- 1), -2.1213203435596433
1: ( 1 -- 2), -2.121320343559643
2: ( 0 -- 2), 1.5
(3) 反力
No.:(節点-方向), 反力
0:( 0 - x), 0.0
1:( 0 - y), 1.5
2:( 2 - y), 1.5(2) 解析モデル２
def datum_exp01()に示す解析モデルと解析結果です。
#出力02
========================================
１．入力データ
(1) 制御データ
ndim, uint_max：[2, 65535]
(2) 共通データ
E：[1.0]
(3) 節点座標
節点数：5
0: [x: 0.00000, y: 0.00000]
1: [x: 100.00000, y: 100.00000]
2: [x: 200.00000, y: 0.00000]
3: [x: 300.00000, y: 100.00000]
4: [x: 400.00000, y: 0.00000]
(4) 部材
部材数：7
0: (node: 0 -- 1), [A: 100.00000]
1: (node: 1 -- 2), [A: 100.00000]
2: (node: 2 -- 3), [A: 100.00000]
3: (node: 3 -- 4), [A: 100.00000]
4: (node: 0 -- 2), [A: 100.00000]
5: (node: 2 -- 4), [A: 100.00000]
6: (node: 1 -- 3), [A: 100.00000]
(5) 境界条件
境界条件を与える節点の数：2
[境界条件：fix:0, free:1]
0: (node: 0), [ 0, 0]
1: (node: 4), [ 1, 0]
(6) 外力
外力が作用する節点の数：2
0: (node: 1), [x: 0.00000, y: -3.00000]
1: (node: 3), [x: 0.00000, y: -3.00000]
========================================
２．計算結果
(1) 節点変位
節点: 変位[x:], [y:], [xy:]
0: [x:0.0], [y:0.0], [xy:0.0]
1: [x:8.999999999999996], [y:-17.485281374238568], [xy:19.665580711901036]
2: [x:5.999999999999996], [y:-20.48528137423857], [xy:21.345883747967093]
3: [x:2.999999999999996], [y:-17.48528137423857], [xy:17.7407740737628]
4: [x:11.999999999999995], [y:0.0], [xy:11.999999999999995]
(2) 部材応力
部材: (節点0--節点1), 部材応力
0: ( 0 -- 1), -4.242640687119284
1: ( 1 -- 2), 2.51214793389404e-15
2: ( 2 -- 3), 0.0
3: ( 3 -- 4), -4.2426406871192865
4: ( 0 -- 2), 2.999999999999998
5: ( 2 -- 4), 2.9999999999999996
6: ( 1 -- 3), -3.0
(3) 反力
No.:(節点-方向), 反力
0:( 0 - x), 1.3322676295501878e-15
1:( 0 - y), 2.999999999999999
2:( 4 - y), 3.0４．改訂履歴
以下、改訂履歴です。2020.3.11にVer.1.03に改訂しました。
改訂履歴の表示・非表示の切り替え
※ブラウザによっては最初から表示されてしまいます。(Google Chrome推奨)
(2019.11.19) Ver.1.01
プロトタイプ
(2020.3.9) Ver.1.02
Ver.1.01：共通
Numpy配列の表記をa[i][j]タイプからa[i, j]タイプに変更
Ver.1.01：44行目
ndimの部分はディメンジョンを表すものではないため2とする
indv = np.zeros((nel, ndim), dtype = 'uint16')
→indv = np.zeros((nel, 2), dtype = 'uint16')
Ver.1.01：194行目
ek_global = np.zeros((4, 4), dtype = 'float64')
→不要
Ver.1.01：223, 224行目
rotation[:] = [[c, s, 0., 0.],
[0., 0., c, s]]
→速度重視
rotation[0, 0] = c
rotation[1, 2] = c
rotation[0, 1] = s
rotation[1, 3] = s
Ver.1.01：249行目
return gk
→削除
Ver.1.01：369行目
print(' 部材: (節点1--節点2), 部材応力')
→print(' 部材: (節点0--節点1), 部材応力')
(2020.3.11) Ver.1.03
Ver.1.02：230～250行目
def set_stiff()のgkを関数内でglobal宣言
→引数でgkを渡し、returnでgkを返す。
Ver.1.02：412行目
set_stiff(iii, ek_global, indv)
→gk = set_stiff(iii, ek_global, indv, gk)
Ver.1.02：281～299行目
def del_stiff_fix()のf_load, gkを関数内でglobal宣言
→引数でf_load, gkを渡し、returnでf_load, gkを返す。
Ver.1.02：412行目
del_stiff_fix(nmat, n_free, ind_free)
→f_load, gk = del_stiff_fix(nmat, n_free, ind_free, f_load, gk)
４．関連記事
以下、平面トラスの有限要素法サンプルコードの解説記事です。
私が実際に購入した教材のご紹介
以下、私が実際に購入したPythonの教材をまとめてみました。 Pythonを学習する上で、少しでもお役に立つことができればうれしいです。
・Python♪私が購入したPythonの書籍のレビュー
・UdemyのPythonの動画講座を書籍を買う感覚で購入してみた
その他
以下、私が光回線を導入した時の記事一覧です。
(1) 2020年「光回線は値段で選ぶ」では後悔する ←宅内工事の状況も説明しています。
(2) NURO光の開通までWiFiルーターを格安レンタルできる
(3) NURO光の屋外工事の状況をご紹介。その日に開通！
(4) 光回線開通！実測するとNURO光はやっぱり速かった
(5) ネット上のNURO光紹介特典は個人情報がもれないの？
