第2章は第1章を受けて簡単な3角形一次要素のプログラムを組むというものである。仕様は次のようになる。
物性値は一つだけ使用可能。
三角形一次要素を使用する
平面ひずみ問題を解く
境界条件は節点荷重のみ使用可能。拘束条件は任意の数値を入力できる。
計算結果(変位・反力・応力・歪)を出力する。
この章の最初の部分にあるべき仕様については全く記載がない。何故書かない。その昔読んだ本で、大学の授業の悪い点として目的を教えないというのがあった。書籍にする以上、目的(ゴール)を書かないでどうする。
実際にプログラムを書く前に仕様を決める必要がある。この章のプログラムなら、先に挙げたものでも十分だろう。まず大まかな仕様を決め、次にロジック(アルゴリズムともいう)を決める必要がある。アルゴリズムと言わないのは、定義された式を如何にプログラムに落とし込むか考える必要があるからである。この部分はそれこそ書く前に考えておく必要がある。考えながら書いてはならない。混乱の元である。
次にプログラムの基礎として構造化プログラムの話やデーター構造の話もない。この分野で最も軽んじられているのがデーター構造で、構造体を使いたくないのか、使えないからなのか、配列で代用している例まで散見される。構造体は初心者には難しいと思われるが、これを超えてくれないとマトモなプログラミングスキルは身につかない。
構造化プログラミングというのは、Dijkstraが提唱したもので、プログラムが次の要素から構成されれていることを言う。
順次実行
プログラムは上から下に処理される。
反復処理
要するにループのことである。ある条件下の元で同じ処理を実行することを言う。
条件分岐
通常のプログラミング言語のif文。ある条件で処理が複数に分かれることを言う。
手続き
本文吹き出しにもあるが、小さなまとまり(サブプログラム)だが、そこには「一つのことを処理する」という意図が必要である。逆に言うとあれもこれもやらないことである。
P28にはこのプログラムのフローチャートなるものが載っている。確か情報関連ではフローチャートはやめましょうと言われていたらしいが、単なる2度手間のために労力を割くのは詰まらないだろう。これくらいなら態々図にしなくても箇条書きで十分だろうし、それにサブルーチン名もいらない。
この後グダグダと説明が続くが、そのようなものは先の箇条書き版で「解説」すれば十分である。文章の構成というかどうしたら理解しやすくできるかという工夫が足らない。
この説明の後にメインルーチンの完成版が登場する。次のとおりである。
' ------------------
メインプロシージャ―
' ------------------
Sub simple_fem() ' メインプロシージャ―
Call initialize ' 配列を初期化
Call make_D ' Dマトリックスの作成
Call make_B ' Bマトリックスの作成
Call make_Ke ' 要素剛性マトリックスの作成
Call make_K ' 全体剛性マトリックスの作成
Call set_boundary_condition ' 境界条件処理
Call solve ' 連立方程式を解く
Call make_reaction ' 節点反力を計算
Call make_strain_element ' 要素ひずみを計算
Call make_stress_element ' 要素応力を計算
End Sub
一言で言ってコメントがしつこい。 コメントをなるべく入れるようにしようと書いてる本があるが、ただのオウム返しならそんなものは要らない。特にサブルーチンのヘッダーにコメントを入れてはならない。サブルーチンが何をしているのかを説明する場合、サブルーチンの前に入れる。サブルーチンの中に書いてはならない。その部分でプログラムが切れるので、読むのに骨が折れ、結局コメントそのものが無視されるのがオチである。
またサブルーチンの名称が中途半端だったり、余計なものがついていたりして不適切である。Initializeは配列の初期化とあるが、それならinit_arrayとした方がいいだろう。この後行列演算用のサブルーチンを呼んでいるが、この4つのルーチンの名前には_matを付けたほうがいい。意地悪な人がmake_Dだけ見たら、これは変形速度テンソルの計算をするんですかと聞いてくるだろう。同じことはmake_Bにも言える。正確を期するとここはmake_Iso_PSE_De_matとでもした方がいいだろうが、この例題ではそこまでする必要はないだろう。
次のmake_Kもだが、make_Keと名前が似通って紛らわしい。このルーチンは名前をmake_SysK_matとでもした方がいいだろう。これであれば要素剛性の計算を意味するmake_Ke_matと間違える心配はない。最もこの部分の問題は別のところにあるのだが。
次のset_boundary_condition以下を見るとなんというか、ポリシーの欠如が見て取れる。どうせその前のルーチンの呼び出しで略号を使っているのだから、set_BCで十分だろう。solveはこのままでいいとして、make_reactionは何を計算するのだろうか。コメントを見ると反力を計算するとある。また、好みの問題かもしれないが、makeは作るであって計算ではない。行列は作成だからmakeでいいかも知れないが、計算するならばcalcrationを用いるべきだろう。むろんこのままでは長いから、calcに短縮する。
以上をまとめると、make_reactionではなく、calc_reaction_forceとするべきだ。片方はフルスペルで書いてもう片方は省略するというのは一貫性の欠如と言われても反論できまい。
最後に呼び出される2つの要素内部の量についても同じことがいえる。FEMでは基本的に積分点であの値を出力するからelementは不要である。またmake_reactionと同じ理由からmakeではなくcalcを用いる。
以上をまとめるとメインルーチンは次のようになる。calc_reaction_forceはさらに名前を短くしている。
' 単純なFEMプログラムの実現
Sub simple_fem()
Call init_array
Call make_D_mat
Call make_B_mat
Call make_Ke_mat
Call make_SysK_mat
Call set_BC
Call solve
Call calc_react_force
Call calc_strain
Call calc_stress
End Sub
名前の省略について説明しておく、名前は基本的に英語ベースで書くが長い単語を書くというのはなかなか面倒なのでどこかで後ろを切るということがある。私の場合、以下の指針を採用している。
発言してみて座りのいいところで切る。
4文字程度にする
単語の長さが5~6文字程度ならそのまま採用する。
各単語は下線でつなぐ。
なるべく元の変数名を反映するようにする。
基本がこれであるが、例外はどこにでも存在する。三角形の場合triangleだからtriで切るより他ない。efem.cではtriaにしているが、これは構造体の初期化の時、四角形要素の方と字数を合わせるためにこのようにしている。
次がデーター構造である。単純に言うと関連するデーターを組織立って使うことで、FEMであれば、例えば、節点データーを定義する際、節点番号と座標値、必要なら局所座標系のインデックスをひとまとめにすることを言う。少なくとも前2者はどこのソルバーでも共通して入力するようになっている。ならば、プログラムでもこれらをひとまとめにしたら使い勝手がいいのではとなるわけだが、ソルバーのみだと実はあんまり有難味はない。
次に要素だが、これも構造体にした方が有難味が大きいのではないかということになる。旧FORTRANのようにipropsなどという惨めな代用品を使うよりデーター名でアクセス出来たほうが余程スマートだし、可読性も向上する。
因みにfem.cでの要素の定義は次のようにしている。
/* 要素データー */
typedef struct {
int number;
PartData *part;
ElementInfo *info;
int *conn;
ElementCalcInfo *cinfo;
} Element;
この構造体は上から順に要素番号、パーツデーターへのポインタ、要素情報へのポインタ、コネクティビティ(実際は整数配列)、最後に計算中の情報を保持するためのデーターへのポインタ(実際は配列)が収められている。作成する基準が違うからデーターの定義の仕方も異なるが、実際のプロダクト並みの能力を持たせ、種々の要素を統一的に使用するという目的からこのようなデーターを組んでいる。
実際に構造体を使う理由は処理速度の向上ではなく、データーを一つのパッケージとして扱うためである。こうすれば個別に変数を定義する必要はなく、構造体そのものをルーチンに渡せば済む。つまり見かけの嵩が減るわけである。
だが、このプログラムではそんなことはしていない。そんな概念すらないようである。VBAには構造体はもとよりクラスまであるというのにである。
このプログラムのトップレベルで宣言されている変数や定数の宣言部分を掲載する。
Option Base 1 ' 配列のインデックスを1から始める
Option Explicit ' 変数の定義を強制する
' ----
' 定数
' ----
Const NODES As Integer = 6 ' 全節点数
Const ELEMENTS As Integer = 4 ' 全要素数
Const COMPONENTS As Integer = 3 ' ひずみと応力の成分数
Const NODES_TRIA3 As Integer = 3 ' 三角形1次要素の節点数
Const DOF_NODE As Integer = 2 ' 節点自由度
Const DOF_TOTAL As Integer = DOF_NODE * NODES ' 全体自由度
Const DOF_TRIA3 As Integer = DOF_NODE * NODES_TRIA3 ' 要素自由度
Const THICKNESS As Double = 1# ' 要素の厚さ
Const YOUNG As Double = 210000# ' ヤング率
Const POISSON As Double = 0.3 ' ポアソン比
' --------------------
' モジュールレベル変数
' --------------------
Dim x(NODES) As Double ' 節点のx座標配列
Dim y(NODES) As Double ' 節点のy座標配列
Dim connectivity(ELEMENTS, NODES_TRIA3) As Integer ' 要素内節点配列
Dim D(COMPONENTS, COMPONENTS) As Double ' Dマトリックス
Dim B(ELEMENTS, COMPONENTS, DOF_TRIA3) As Double ' Bマトリックス
Dim area_element(ELEMENTS) As Double ' 要素面積の配列
Dim Ke(ELEMENTS, DOF_TRIA3, DOF_TRIA3) As Double ' 要素剛性マトリックス
Dim K(DOF_TOTAL, DOF_TOTAL) As Double ' 全体剛性マトリックス
Dim U(DOF_TOTAL) As Double ' 全体変位ベクトル
Dim F(DOF_TOTAL) As Double ' 全体荷重ベクトル
Dim Um(DOF_TOTAL) As Boolean ' 拘束目印ベクトル
Dim Kc(DOF_TOTAL, DOF_TOTAL) As Double ' 計算用全体剛性マトリックス
Dim Fr(DOF_TOTAL) As Double ' 節点反力ベクトル
Dim strain_element(ELEMENTS, COMPONENTS) As Double ' 要素ひずみベクトル
Dim stress_element(ELEMENTS, COMPONENTS) As Double ' 要素応力ベクトル
なんというか微妙な感が漂う部分である。まず全節点数や全要素数を意味する定数NODESやELEMENTSはNTNODE(Number of Total NODEs)とNTELEM(Number of Total ELEMents)とした方がいいだろう。先頭にNを付けた変数は個数を表すルールを一様に適用すればいい。次のCOMPONENTSはNCOMPでいいだろう。それに抑々固有名詞の類は型名に使うべきである。
DOF_NODEはMAXDOFとする。この項は節点自由度ではなく、最大自由度とした方がいい。こうしておけば、自由度の方向は置くにしてもこれが自由度の個数の上限ということが一目でわかる筈である。
DOF_TOTALは初めの例にならってNTDOFでいいだろう。問題はDOF_TRIA3である。この定数は3角形1次要素の節点数を意味する変数だが、これだけで完結し、その後要素の種類を追加する気が無いなら、NENODE(Number of NODE at Element)で十分である。
複数の要素形状(とどのつまり三角形と四角形)を使えるようにするにはどうしたらいいかという問題があるが、今のところ本題ではないから置くとする。一つだけ指摘しておくと、この要素は後で出てくるアイソパラメトリック要素と相性が悪いのである。そのため実用レベルを目指すなら三角形1次要素もアイソパラメトリック定式化を行った方がよい。
次が変数の宣言であるが、物性値は1種類だからこれでいいとして、座標値は一つの配列とするべきである。この様な要素では大した問題にはならないが、アイソパラメトリック要素を実装するときに苦労する。また、基本的に大域レベルで1～2文字の変数を使うものではない。探すのに苦労する。最悪でも4文字からとするべきである。
Umも問題があるがとりあえず安直と言える名前を付けるのはやめにした方がいい。たとえば、DOF_FLAGSの方がよかっただろうに。同じことはKcにも言える。
最後のstrain_elementとstress_elementの_elementは自明だから消去して良い。FEMの計算では通常要素やひずみの値は積分点で取る。この要素は要素全体で応力・ひずみが一定である。いずれにしても要素の値を求めているから態々elementをつける意味がない。
だがこのメインルーチンと変数の設計にはもっと大きな問題がある。まず、Bマトリックスを丸ごと保存しているが、データー構造を工夫すれば三角形一次要素であれば6個に減らせるし、ひずみの計算で再使用する際も簡単な計算で行うことが出来る。それより大きな問題は要素剛性マトリックスで、実はこの配列は一つ分だけあれば十分である。何故かというと、このプログラムではBマトリックス・要素剛性マトリックス・全体剛性マトリックスを順次計算するが、実用的な観点から言うと、全体剛性マトリックスの計算を行う中でBマトリックス・要素剛性マトリックスの計算を順次行う。解かりやすいようにというが、あまりにも馬鹿にしすぎやしないか。FEMのプログラムを組む方に指摘しておくが、「初めての有限要素法プログラム」なら致し方ないが、この後もこのような方法を墨守すると周囲から笑われる程度のものであると断っておく。
この部分の実例を以下に掲げる。
typedef void (*FP_calc_B) (double *, const Element *, const double *, const double *, int);
void calc_K_mat (FILE *fout)
{
int i, j;
int nvstat = etinfo [prob].nvstat;
FP_calc_B calc_B [] = { NULL, calc_PL_B_mat, calc_PL_B_mat, calc_AXSOL_B_mat, NULL, };
create_K_info (fout);
calc_load ();
for (i = 0; i < ntelem; i++) {
Element *e = elem + i;
ElemIntegralInfo *info = e->info->info1;
double *N = info->N;
double *H = info->dN;
double *D = e->part->init_mat;
int ndim = etinfo [prob].ndim, nnode = info->nnode;
for (j = 0; j < info->ntint; j++) {
(*calc_B [prob]) (__B__, e, N, H, j);
attach_B_mat (e, D, j);
N += nnode;
H += nnode * ndim;
}
}
dump_pline (OFF, OFF);
}
この関数(Cではサブルーチンも関数と称する)では二重のループから成るが、外周のループは全要素についてのもので、内周のループが全積分点に関するものである。(*calc_B [prob]) (__B__, e, N, H, j)が要素のひずみ変位マトリックスの作成ルーチンの呼び出しで、次のattach_B_mat (e, D, j)が全体剛性マトリックスに要素剛性マトリックスの繰り込みを行う関数である。要素剛性マトリックスはこの関数で作られることになるが、実はこの中では要素剛性マトリックス全体を明示的に作っていない。それは置くとして、実際の全体剛性マトリックスの作成はこのようなものであり、一つ一つ実行するというようなものではない。
もし効率化というなら、アルゴリズムの時点でどうやったら早くなるか検討してから実装したほうがよいが、それ以前の問題としてあからさまに非効率かつ幼稚とすら言える方法を紹介するのは感心しないのである。
訂正：20190907 文章を一部変更