有限要素法・アルゴリズム - 超音波流体屋のプログラム備忘録

六面体要素の例でアルゴリズム処理の並びを整理する。
プログラム上での便宜のために導出で用いなかった記号を用いる。

定数行列の準備

${\cal C}_{ijkl} = \lambda \delta_{ij}\delta_{kl} + \mu (\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}) $
$ S_{kn} = \left[ \begin{array}{cccccccc} 1 & -1 & -1 & 1& 1& -1& -1& 1 \\ 1 & 1 & -1 & -1& 1& 1& -1& -1 \\ 1 & 1 & 1 & 1& -1& -1& -1& -1 \end{array} \right]$

メッシュ

以下のように作成済みとする。
$ x_{ip} = \left[ \begin{array}{ccc} x_1 & & x_p\\ y_1 & \cdots & y_p \\ z_1 & & z_p \end{array} \right] ,\ e_{nq} = \left[ \begin{array}{ccc} p_{11} & & p_{1q}\\ \vdots & \cdots & \vdots \\ p_{n1} & & p_{nq} \end{array} \right] $

自由度・荷重

拘束・外力込みで定義しておく、このページでは節点における自由度を先にカウントする。
自由度に対する時間微分などあれば同様。
$ u = ^t\left[ \begin{array}{ccc} u_{x1}, u_{y1}, u_{z1}, & \cdots, & u_{zp} \end{array} \right]$
$ f = ^t\left[ \begin{array}{ccc} f_{x1}, f_{y1}, f_{z1}, & \cdots, & f_{zp} \end{array} \right]$
全自由度の配列サイズは$L_D = L_p L_d$ (配列pのサイズを $L_p$ = length(p) とする)。

要素内M,K行列の作成

結構重たい処理のため、並列化ができる今の御時世最低smpで同時処理するのが望ましい。
parfor: は並列化推奨ループを示す。

parfor q = 1..q(各要素に対するループ)

• $N^q(n, s), D^{q\prime}(d,n,s), G(d,d,s), detJ(s)$: メモリの動的確保

parfor s = 1..n(数値積分各要素に関するループ)

1. $\xi_{ks} = S_{ks}/\sqrt{3} $　　
2. $N_n = \prod_k \frac{1 + S_{kn} \xi_{ks}}{2} $
3. $D^\prime_{jn} = \frac{\partial N_n}{\partial \xi_j} = S_{nj} \prod_k \frac{1 + \bar{\delta}_{jk} S_{nk} \xi_{ks}}{2}$
4. $J_{ij} = x_{in} D^\prime_{jn} $ 　　[loop: n]
5. $G_{ij} = J_{ij}^{-1} $
6. $dV = det J_{ij} $

end(index s)

end(要素 q)

並列化パラメータの関係で一旦要素ループをここで切る。

parfor q = 1..q(各要素に対するループ)

• $M^q(n,n), K^q(n,n), {\cal M}^q(d,d,n,n), {\cal K}^q(d,d,n,n)$行列のメモリの動的確保・0フィル

parfor m = 1..Σn, n = 1..Σn(K行列等各要素に対するループ)

• $D^{n,m}(d), {\cal B}^{n,m}(d,d,d)$行列のメモリ確保
• for s = 1..n(数値積分): リダクション(1メモリへの複数回アクセス)が噛むので並列化しないのが安泰

1. $D_{im}, D_{in} = N_{n,i} = D^\prime_{jn} G_{ij} $　　[loop: j]
2. ${\cal B}_{uvjm}, {\cal B}_{klin} = \frac{1}{2} \left[ \delta_{ik} D_{ln} + \delta_{il} D_{kn} \right]$
3. $K^q_{mn} += D_{im} D_{in} dV$　　[loop: i]
4. ${\cal K}^q_{ijmn} += {\cal B}_{uvjm} {\cal C}_{uvkl} {\cal B}_{klin} dV $ [loop: k,l,u,v ]
5. $M^q_{mn} += N_m N_n dV$
6. (意図してN行列ベースの質量行列を作成したければ)${\cal M}^q_{ijmn} += \rho \delta_{ij} N_m N_n dV$

• end (index s)

end(index m, n)

end(要素 q)

• $N^q(n, s), D^{q\prime}(d,n,s), G(d,d,s), detJ(s)$: は以後使わないので解放する。

拘束・非ゼロ情報の登録

陰解法など全体行列が必要な場合には、合間に(CPUなど)単プロセスで自由度IDに関する整備をしておく。
メッシュ登録の段階くらいで、前もって隣接節点状況を登録しておくこと。
free = int[L_D], free_count = 0
メモリ節約のため、密自由度座標を疎行列インデックスへの変換を記録するidxは帯行列とする。
$idx = \textrm{int}\left\{ L_D \times L_d\max[L_{p,a}] \right\}$, idx_count = 0;

for m = 1..Lp(各節点に対するループ)

for i = 1..Ld(節点mの自由度)

• $r = i + L_d m$
• if free(m,i) then

　　free(r) = free_count
　　free_count++

• else

　　free(r) = -1

• end
• for a = 1..max L(p,a)(隣接節点に対するループ, $a > L(p,a)$ ならリターン )
• for j = 1..Ld(節点aの自由度)
• $n = m.adjacent(a)$ : 隣接節点の番号
• $c = j + L_d a$
• if free(m,i) and free(n,j) then

　　idx(r,c) = idx_count
　　idx_count++

• else

　　idx(r,c) = -1

• end
• end(index j)
• end(隣接節点 a)

end(index i)

end(節点 m)

全体疎行列

拘束自由度を除いた全体行列を作成する。
メモリの確保は $(K,M)$ = double [idx_count] $(row, col)$ =int [idx_count]
及び $\tilde{f} =$ double [free_count]

parfor m = 1..Lp(各節点に対するループ)

parfor i = 1..Ld(節点mの自由度)

1. $r = i + L_d m$
2. $l = free(r)$
3. $\tilde{f}(idx_f) = f(r) \hspace{2em} if\ l \neq -1$

• for a = 1..max L(p,a)(隣接節点に対するループ, $a > L(p,a)$ ならリターン )
• for j = 1..Ld(節点aの自由度)

1. $n = m.adjacent(a)$ : 隣接節点の番号
2. $qlist =m.elem(a)$ : 共有要素リスト
3. $c = j + L_d a$
4. $k = idx(r,c)$
5. $row(k) = r \hspace{2em} if\ k \neq -1$
6. $col(k) = j + L_d n \hspace{2em} if\ k \neq -1 $

• for q = 1..qlist(各要素に対するループ)

1. $ K(k) += {\cal K}^q_{ijmn} \hspace{2em} if\ k \neq -1$
2. 必要なら $M(k) += {\cal M}^q_{ijmn} \hspace{2em} if\ k \neq -1 $
3. $\tilde{f}(l) -= {\cal K}^q_{ijmn} u( j + L_d n) \hspace{2em} if\ l \neq -1 $

• end(要素 q)
• end(index j)
• end(隣接節点 a)

end(index i)

end(節点 m)

最終的な連立方程式に

$ A = K - \omega^2 M $など

連立方程式を解く

$ \tilde{\bf u} = A^{-1} \tilde{\bf f} $
直接法(dgesv_やmatlabならバックスラッシュ)・反復法(CG法他)など適切なルーチンで。
$free$を利用して$\tilde{\bf u}$を配り直して、出力しておしまい。