有限要素法・付録 - 超音波流体屋のプログラム備忘録

有限要素法・導出の付録部分。

積分形

今回の導出で使わないと思うが、積分形ならこんな形をしてる。
(弾): $ \iint {\bf \sigma} d{\bf S} = \iiint \rho (\ddot{\bf u} - {\bf f}) dV $
(ポ): $ \iint \nabla \varphi d{\bf S} = - \iiint \varrho/\varepsilon dV $
(波): $ \iint \nabla \phi d{\bf S} = \iiint (\ddot \phi/c^2 + q) dV $

総和規約

このページでは、テンソル演算において総和記号$\Sigma$は使わない。
$\delta{\bf \epsilon} : {\bf \sigma} = \epsilon_{ij} \sigma_{ij} = \sum_{ij} \epsilon_{ij} \sigma_{ij}$

全体行列の編み方・ボツ

並列計算性が著しく悪い。

• 全自由度: $L_D = L_p L_d$ (配列pのサイズを $L_p$ = length(p) とする)
• 普通は巨大行列になるので、$L_D \times L_D$の疎行列として定義する
• $r, c, M, K, h$: 一次元(空)配列
• $map$ (key:密行列ID)→(value:疎行列ID) : 連想配列 std::map など
• $free$ (key:自由度ID)→(value:非拘束自由度ID) : 連想配列
• $\ll$ = push_back関数: リストに追加を示す

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

for m = 1..Ln(要素内節点に対するループ)

for i = 1..Ld(1節点に対する自由度に対するループ)

• for n = 1..Ln(相手要素内節点に対するループ)
• for j = 1..Ld(相手1節点に対する自由度に対するループ)

1. $p^m = e_{nq}$
2. $p^n = e_{nq}$
3. $p^r = i + L_d p^m$
4. $p^c = j + L_d p^n$
5. $key = p^r + L_D p^c $
6. $idx = map(key) $

• if (idx == NULL) : これまで登録したこと無い要素なら

1. $ map(key) = size(map)$
2. $ r \ll p_r$
3. $ c \ll p_c$
4. $ K \ll {\cal K}^q_{ijmn},\ M \ll {\cal M}^q_{ijmn} $

• else

1. $ K(idx) += {\cal K}^q_{ijmn},\ M(idx) += {\cal M}^q_{ijmn} $

• end
• end(自由度 j)
• end(節点 n)

• if 非拘束自由度(m,i) then $ free(p_r) = size(freemap)$

end(自由度 i)

end(節点 m)

end(要素 q)

拘束付き全体行列・ボツ？

拘束点を間引く必要があるので、一辺$size(free)$のスパース行列として定義する
$r^*, c^*, M^*, K^*, f^*$: 一次元(空)配列

for k = 1..size(r)(全体疎行列に対するループ)

1. $ r^f = free( r(k) ) $
2. $ c^f = free( c(k) ) $
3. next if $r_f = \emptyset \ or \ c_f = \emptyset$ 拘束されているなら登録しない
4. $ r^* \ll r_f $
5. $ c^* \ll c_f $
6. $ K^* \ll K(k),\ M^* \ll M(k)$ など

end(index k)

for k = 1..L_D(自由度に対するループ)

1. $ r^f = free(k) $
2. $ f^* \ll \tilde{f}(k) $ if $r_f \neq \emptyset $

end(index k)

最終的な連立方程式に

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