前回に引き続いて2.2節の批評に移る。というか前回書いているときに気が付かなかったのだが、前回の大域変数の定義は2.2節「変数と配列」の部分に当たる。
この部分の文章の組み立てとしてはとりあえず宣言する大域変数の説明を行い、それからモデルの説明をするのが筋というものだろう。
さて、このプログラムのコンセプトとして、データーをプログラムにハードコードして、お手軽に経験してもらうというのがある。この是非は置くとしても、やり方がまずい。文中のプログラム2.4を掲載する。
Sub initialize()
' 要素内節点順を設定
connectivity(1, 1) = 1: connectivity(1, 2) = 2: connectivity(1, 3) = 5
connectivity(2, 1) = 2: connectivity(2, 2) = 3: connectivity(2, 3) = 4
connectivity(3, 1) = 2: connectivity(3, 2) = 4: connectivity(3, 3) = 5
connectivity(4, 1) = 1: connectivity(4, 2) = 5: connectivity(4, 3) = 6
' 節点座標配列の設定
x(1) = 0#: y(1) = 0#
x(2) = 1#: y(2) = 0#
x(3) = 2#: y(3) = 0#
x(4) = 2#: y(4) = 1#
x(5) = 1#: y(5) = 1#
x(6) = 0#: y(6) = 1#
End Sub
何が何だか分からない、というより理解したくない。この方法の欠点は値を直接配列に設定しているため、プログラムをよく読まないといけないことなのだが、これで間違いがあったら目も当てられない。
このようなルーチンを作るくらいなら手間をかけて順にコネクティビティ(要素における節点の並びのこと)をセットするset_connや座標値を設定するset_coordといった手続きを作ってそこで設定したほうがいい。
例えばこのようになる。
Sub set_coord(byval node as integer, byval xc as double, byval yc as double)
x(node)=xc: y(node)=yc
End Sub
sub set_conn(byval elem as integer, ParamArray conn() As Variant)
dim i as integer
for i=1 to ubound(conn)
connectivity(elem,i)=conn(i)
next i
end sub
set_coordは節点番号と座標値を渡し、それをx,yに渡す。set_connは初心者向きではないが、最後のconnが可変長引数となっており、これを使ってコネクティビティのセットをする。このサブルーチンのuboundは配列の上限値を返す。
このようなサブルーチンを頑張って作れば、あとは機械的にこのルーチンをを呼ぶだけで値の設定が終了する。
そしてInitializeは次のようになる。あの意味不明なルーチンはこのように単純になる。サブルーチンを追加しただけだが、このような手間をかけてプログラムの可読性を上げたほうがいいだろう。
Sub initialize()
' 要素内節点順を設定
set_conn 1, 1, 2, 5
set_conn 2, 2, 3, 4
set_conn 3, 2, 4, 5
set_conn 4, 1, 5, 6
' 節点座標配列の設定
set_coord 1, 0, 0
set_coord 2, 1, 0
set_coord 3, 2, 0
set_coord 4, 2, 1
set_coord 5, 1, 1
set_coord 6, 0, 1
End Sub
次が2.3節「要素剛性マトリックス」で例の誤った定式化が載っている以外前置きについてはこれでいいだろう。
2．3.1節は平面ひずみのDマトリックスを作成するmake_Dの説明である。
1: Sub make_D()
2: Dim coef As Double ' ?}?g???b?N?X??????????W??
3: coef = YOUNG / (1# - 2# * POISSON) / (1# + POISSON)
4:
5: D(1, 1) = coef * (1# - POISSON) ' D?}?g???b?N?X???
6: D(1, 2) = coef * POISSON
7: D(1, 3) = 0#
8: D(2, 1) = D(1, 2)
9: D(2, 2) = coef * (1# - POISSON)
10: D(2, 3) = 0#
11: D(3, 1) = D(1, 3)
12: D(3, 2) = D(2, 3)
13: D(3, 3) = coef * (1# - 2# * POISSON) / 2#
14: End Sub
文章の方は言いがプログラムに無駄な代入があり、それに追い打ちをかけるようにつまらない吹き出しがある。本書の吹き出しにて曰く、
プログラムの8,11,12行目を個別に記述しても構いませんが、既に計算された変数を再利用をするとプログラムミスが減少し、保守も容易になります。
保守性云々はその通りであるが、「違う、そうじゃない」。
プログラミング上達のコツの一つに同じことを書くなというのがある。このことを言うのならば定数で済む部分ではなく計算を行っている式だけで説明した方がよかっただろう。実際に再利用するべきは5行目のD(1,1)と6行目のD(1,2)であり、この結果をD(2,2)とD(2,1)にそれぞれ代入する。成分が0のところはそのまま値を代入する。というか、態々変数を代入するというのは無駄ですらある。これなら複文にしても問題はないから一行にまとめておく。次に一時変数のcoefだが、一文字で迷いようがないからfでいいだろう。
以上をまとめると次のようになる。
1: Sub make_D()
2: Dim f As Double
3:
4: f = YOUNG / (1 - 2 * POISSON) / (1 + POISSON)
5: D(1, 1) = f * (1 - POISSON)
6: D(1, 2) = f * POISSON
7: D(2, 1) = D(1,2)
8: D(2, 2) = D(1,1)
9: D(1, 3) = 0: D(2, 3) = 0: D(3, 1) = 0: D(3, 2) = 0
10: D(3, 3) = YOUNG * (1 + POISSON) / 2
11: End Sub
この後の節で、吹き出しの中で「Dマトリックスが一般的に非対称になる」というのがあるがこれは不要である。本来の定式化を行えばDは転置にならないからそのまま計算することになるからである。それにこの部分は初心者向けとは思えない。
素朴な疑問だが構造解析分野でDマトリックスが非対称になるという問題はどれくらいあるのだろうか。彼らによると一般的にな、とあるから相当あるのだろう。筆者は過分にして知らないが、非線形(弾塑性)解析では2ケースほどあるのは知っている。どれくらいあるんでしょうかね。
2.3.2のBマトリックスの作成についてはとりあえずこのままでいいだろうが、Bマトリックスそのものを要素分格納しなくても計算した項目だけ格納すれば必要な領域は1/3に縮小することを指摘しておく。
参考までの原著のBマトリックスの計算ルーチンを掲載する。
Sub make_B()
Dim e As Integer ' 要素カウンター
For e = 1 To ELEMENTS ' 要素ごとに処理
' 節点座標をローカル変数に代入
Dim x1 As Double, y1 As Double
Dim x2 As Double, y2 As Double
Dim x3 As Double, y3 As Double
x1 = x(connectivity(e, 1)): y1 = y(connectivity(e, 1))
x2 = x(connectivity(e, 2)): y2 = y(connectivity(e, 2))
x3 = x(connectivity(e, 3)): y3 = y(connectivity(e, 3))
' 要素面積を計算
area_element(e) = (x1 * y2 - x1 * y3 + x2 * y3 - x2 * y1 + x3 * y1 - x3 * y2) / 2
Dim coef As Double ' マトリックス成分に共通な係数
coef = 1# / area_element(e) / 2#
B(e, 1, 1) = coef * (y2 - y3) ' Bマトリックスを作成
B(e, 1, 2) = 0#
B(e, 1, 3) = coef * (y3 - y1)
B(e, 1, 4) = 0#
B(e, 1, 5) = coef * (y1 - y2)
B(e, 1, 6) = 0#
B(e, 2, 1) = 0#
B(e, 2, 2) = coef * (x3 - x2)
B(e, 2, 3) = 0#
B(e, 2, 4) = coef * (x1 - x3)
B(e, 2, 5) = 0#
B(e, 2, 6) = coef * (x2 - x1)
B(e, 3, 1) = B(e, 2, 2)
B(e, 3, 2) = B(e, 1, 1)
B(e, 3, 3) = B(e, 2, 4)
B(e, 3, 4) = B(e, 1, 3)
B(e, 3, 5) = B(e, 2, 6)
B(e, 3, 6) = B(e, 1, 5)
Next e
End Sub
2.3.3が要素剛性の計算の説明になる。最も初心者向きなのか馬鹿々々しさを感じる書き方を行っている。プログラム2.14がそれである。取り敢えず掲載するのでじっくり味わっていただきたい。
Sub make_Ke()
Dim Dt(COMPONENTS, COMPONENTS) As Double ' Dの転置マトリックス
Dim Bt(DOF_TRIA3, COMPONENTS) As Double ' Bの転置マトリックス
Dim BtDt(DOF_TRIA3, COMPONENTS) As Double ' Bt×Dtのマトリックス
Dim e As Integer ' 要素番号カウンター
Dim r As Integer, c As Integer ' マトリックスの行と列のカウンター
Dim m As Integer ' マトリックスの掛け算用カウンター
For r = 1 To COMPONENTS ' Dの転置マトリックスDtを作成
For c = 1 To COMPONENTS
Dt(c, r) = D(r, c)
Next c
Next r
For e = 1 To ELEMENTS ' 要素ごとに処理
For r = 1 To COMPONENTS ' Bの転置マトリックスBtを作成
For c = 1 To DOF_TRIA3
Bt(c, r) = B(e, r, c)
Next c
Next r
For r = 1 To DOF_TRIA3 ' Bt×Dt
For c = 1 To COMPONENTS
BtDt(r, c) = 0#
For m = 1 To COMPONENTS
BtDt(r, c) = BtDt(r, c) + Bt(r, m) * Dt(m, c)
Next m
Next c
Next r
For r = 1 To DOF_TRIA3 ' Keの計算
For c = 1 To DOF_TRIA3
Ke(e, r, c) = 0#
For m = 1 To COMPONENTS
Ke(e, r, c) = Ke(e, r, c) + BtDt(r, m) * B(e, m, c)
Next m
Ke(e, r, c) = Ke(e, r, c) * area_element(e) * THICKNESS
Next c
Next r
Next e
End Sub
これを見た時の第一印象は何ぞこれ、だった。本来の定式化では$D^T$の計算は不要であるからこの部分は丸ごとカットできる。そして$B^T$の計算もカットできる。このようなものは一々プログラムを書かずとも、掛け算を実行する際に添え字を入れ替えるだけで同じことが出来る。ループカウンターに一文字の変数を使うのはお約束だが、これではわかりづらいし、コメントも役に立っていない。更に公式に従うべきなのでここも修正する。
この部分を修正するとこのようになる。
' 要素剛性を計算する
Sub make_Ke()
Dim BtDt(DOF_TRIA3, COMPONENTS) As Double ' BtとDの積
Dim i As Integer, j As Integer, k As Integer, As Integer
For i = 1 To ELEMENTS
For j = 1 To DOF_TRIA3
For k = 1 To COMPONENTS
BtDt(j, k) = 0
For l = 1 To COMPONENTS
BtDt(j, k) = BtDt(j, k) + e(i, l, j) * Dt(l, k)
Next l
Next k
Next j
For j = 1 To DOF_TRIA3
For k = 1 To DOF_TRIA3
Ke(e, j, k) = 0
For m = 1 To COMPONENTS
Ke(e, j, k) = Ke(e, j, k) + BtDt(j, l) * B(e, l, k)
Next m
Ke(i, j, k) = Ke(i, j, k) * area_element(i) * THICKNESS
Next k
Next j
Next i
End Sub
行数は42行から25行に減った。アルゴリズムの工夫と、公式に従っただけでこれである。どれだけ非効率か分かるだろう。この分野に顕著だが、長ったらしいだけのサブルーチンを書くのを「美徳」と持っている輩が多いが、現在要求されていることは「よりコンパクトに」である。
しかしこの程度で満足してはいけない。行列の積や転置を求める関数を作成し、それらを呼び出せば、このサブルーチンの内容をもっとすっきりさせることが出来ただろうが、この設計では不可能である。VBAは部分配列が取れない為である。
2日ほどで作成した単純なFEMプログラムからこの部分を抜粋する。
'行列の積を求む
Function matmul(a() As Double, B() As Double)
Dim res() As Double
Dim i As Integer, j As Integer, k As Integer
ReDim res(UBound(a, 1), UBound(B, 2))
For i = 1 To UBound(a, 1)
For j = 1 To UBound(B, 2)
res(i, j) = 0
For k = 1 To UBound(B, 1)
res(i, j) = res(i, j) + a(i, k) * B(k, j)
Next k
Next j
Next i
matmul = res
End Function
'転置行列を求む
Function transpose(a() As Double)
Dim res() As Double
Dim i As Integer, j As Integer, k As Integer
ReDim res(UBound(a, 2), UBound(a, 1))
For i = 1 To UBound(a, 2)
For j = 1 To UBound(a, 1)
res(i, j) = a(j, i)
Next j
Next i
transpose = res
End Function
これらの関数はそれぞれ行列の積と転置行列を求める。配列の大きさは先ほど出てきたuboundを用いている。この関数はともにvariant型の変数として定義されている。こうしないとVBAでは配列を返してくれそうにないからこうしている。
このような汎用的な関数やサブルーチンを作成すると、プログラムをいちいち作成する必要がなくなり保守性も向上する。
この関数の使用例は次のとおりである。戻り値の関係上、結果を一旦ローカル変数で受けているが、この時、配列の大きさを指定していないのがミソである。
Sub make_sysK_mat()
Dim i As Integer, j As Integer, k As Integer
Dim B_mat(3, 6) As Double, Ke() As Double, r1() As Double, r2() As Double
Dim nx As Integer, ny As Integer, dx As Integer, dy As Integer
ReDim area(ntelem)
ReDim A_mat(ntelem, 2, 3)
ReDim sysK(ntnode * 2, ntnode * 2)
For i = 1 To ntelem
calc_A_mat i
make_B_mat B_mat, i
r1 = transpose(B_mat)
r2 = matmul(D_mat, B_mat)
Ke = matmul(r1, r2)
For j = 1 To 3
nx = conn(i, j)
For k = 1 To 3
ny = conn(i, k)
sysK(nx * 2 - 1, ny * 2 - 1) = sysK(nx * 2 - 1, ny * 2 - 1) + Ke(j * 2 - 1, k * 2 - 1)*area(i)
sysK(nx * 2 + 0, ny * 2 - 1) = sysK(nx * 2 + 0, ny * 2 - 1) + Ke(j * 2 + 0, k * 2 - 1)*area(i)
sysK(nx * 2 - 1, ny * 2 + 0) = sysK(nx * 2 - 1, ny * 2 + 0) + Ke(j * 2 - 1, k * 2 + 0)*area(i)
sysK(nx * 2 + 0, ny * 2 + 0) = sysK(nx * 2 + 0, ny * 2 + 0) + Ke(j * 2 + 0, k * 2 + 0)*area(i)
Next k
Next j
Next i
End Sub
'要素のAマトリックスを求む
Sub calc_A_mat(i As Integer)
'ヤコビアン(面積)を求む
area(i) = (coord(conn(i,1)*2-1)-coord(conn(i,3)*2-1)) * (coord(conn(i,2)*2)-coord(conn(i,3)*2)) _
- (coord(conn(i,2)*2-1)-coord(conn(i,3)*2-1)) * (coord(conn(i,1)*2)-coord(conn(i,3)*2))
'Aマトリックスの計算
A_mat(i,1,1) = (coord(conn(i,2)*2) - coord(conn(i,3)*2))/area(i)
A_mat(i,1,2) = (coord(conn(i,3)*2) - coord(conn(i,1)*2))/area(i)
A_mat(i,1,3) = (coord(conn(i,1)*2) - coord(conn(i,2)*2))/area(i)
A_mat(i,2,1) = (coord(conn(i,3)*2-1) - coord(conn(i,2)*2-1))/area(i)
A_mat(i,2,2) = (coord(conn(i,1)*2-1) - coord(conn(i,3)*2-1))/area(i)
A_mat(i,2,3) = (coord(conn(i,2)*2-1) - coord(conn(i,1)*2-1))/area(i)
area(i) = area(i) / 2
End Sub
'Bマトリックスを求む
Sub make_B_mat(B() As Double, i As Integer)
Dim j As Integer
For j = 1 To 3
B(1, j * 2 - 1) = A_mat(i, 1, j)
B(1, j * 2 + 0) = 0
B(2, j * 2 - 1) = 0
B(2, j * 2 + 0) = A_mat(i, 2, j)
B(3, j * 2 - 1) = A_mat(i, 2, j)
B(3, j * 2 + 0) = A_mat(i, 1, j)
Next j
End Sub
設計思想からして異なるから、単純な比較はできないものの、いかに効率よくデーターを使うか迄考える必要がある。幸いにして線形解析を行うプログラムは構造が単純なので、このような工夫が生きる余地が大いにある。
2019.9.15 calc_A_matの修正