「有限要素法・アルゴリズム」の編集履歴（バックアップ）一覧はこちら

「有限要素法・アルゴリズム」(2015/07/02 (木) 18:29:06) の最新版変更点

追加された行は緑色になります。

削除された行は赤色になります。

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

*定数行列の準備
$${\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] $$
//添字の連続を避けるため、ある要素qにおけるn番節点のi座標座標を指定するときは $$ x_i(e_{nq}) $$と呼ぶものとする。

*自由度・荷重
拘束・外力込みで定義しておく、このページでは節点における自由度を先にカウントする。
自由度に対する時間微分などあれば同様。
$$ 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: は並列化推奨ループを示す。

**&bold(){parfor q = 1..q}(各要素に対するループ)
-$$N^q(n, s), D^{q\prime}(d,n,s), G(d,d,s), detJ(s)$$: メモリの動的確保
**&bold(){parfor s = 1..n}(数値積分各要素に関するループ)
+$$\xi_{ks} = S_{ks}/\sqrt{3} $$　　
+$$N_n = \prod_k \frac{1 + S_{kn} \xi_{ks}}{2} $$
+$$D^\prime_{jn} = \frac{\partial N_n}{\partial \xi_j} =  S_{nj} \prod_k \frac{1 + \bar{\delta}_{jk} S_{nk} \xi_{ks}}{2}$$
+$$J_{ij} = x_{in} D^\prime_{jn} $$ 　　[loop: n]
+$$G_{ij} = J_{ij}^{-1} $$
+$$dV = det J_{ij} $$
**&bold(){end }(index s)
**&bold(){end }(要素 q)

並列化パラメータの関係で一旦要素ループをここで切る。
**&bold(){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フィル
**&bold(){parfor m = 1..Σn, n = 1..Σn}(K行列等各要素に対するループ)
-$$D^{n,m}(d), {\cal B}^{n,m}(d,d,d)$$行列のメモリ確保
-&bold(){for s = 1..n}(数値積分): リダクション(1メモリへの複数回アクセス)が噛むので並列化しないのが安泰
+$$D_{im}, D_{in} = N_{n,i} =  D^\prime_{jn}  G_{ij} $$　　[loop: j]
+$${\cal B}_{uvjm}, {\cal B}_{klin} = \frac{1}{2} \left[ \delta_{ik} D_{ln}  + \delta_{il} D_{kn} \right]$$
+$$K^q_{mn} += D_{im} D_{in} dV$$　　[loop: i]
+$${\cal K}^q_{ijmn} += {\cal B}_{uvjm}  {\cal C}_{uvkl} {\cal B}_{klin}  dV $$ [loop: k,l,u,v ] 
+$$M^q_{mn} += N_m N_n dV$$
+(意図してN行列ベースの質量行列を作成したければ)$${\cal M}^q_{ijmn} += \rho \delta_{ij}  N_m N_n dV$$
-&bold(){end } (index s)
**&bold(){end }(index m, n)
**&bold(){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;
**&bold(){for m = 1..Lp}(各節点に対するループ)
**&bold(){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
-&bold(){for a = 1..max L(p,a)}(隣接節点に対するループ, $$a > L(p,a)$$ ならリターン )
-&bold(){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
-&bold(){end }(index j)
-&bold(){end }(隣接節点 a)
**&bold(){end }(index i)
**&bold(){end }(節点 m)


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

**&bold(){parfor m = 1..Lp}(各節点に対するループ)
**&bold(){parfor i = 1..Ld}(節点mの自由度)
+$$r = i + L_d m$$
+$$l = free(r)$$
+$$\tilde{f}(idx_f) = f(r)   \hspace{2em} if\ l \neq -1$$
-&bold(){for a = 1..max L(p,a)}(隣接節点に対するループ, $$a > L(p,a)$$ ならリターン )
-&bold(){for j = 1..Ld}(節点aの自由度)
+$$n = m.adjacent(a)$$ : 隣接節点の番号
+$$qlist =m.elem(a)$$ : 共有要素リスト
+$$c = j + L_d a$$
+$$k = idx(r,c)$$
+$$row(k) = r   \hspace{2em} if\ k \neq -1$$
+$$col(k) = j + L_d n   \hspace{2em} if\ k \neq -1 $$
-&bold(){for q = 1..qlist}(各要素に対するループ)
+$$ K(k) += {\cal K}^q_{ijmn}   \hspace{2em} if\ k \neq -1$$ 
+必要なら $$M(k) += {\cal M}^q_{ijmn}  \hspace{2em} if\ k \neq -1 $$
+$$\tilde{f}(l) -=  {\cal K}^q_{ijmn} u( j + L_d n)  \hspace{2em} if\ l \neq -1 $$
-&bold(){end }(要素 q)
-&bold(){end }(index j)
-&bold(){end }(隣接節点 a)
**&bold(){end }(index i)
**&bold(){end }(節点 m)


**最終的な連立方程式に
$$ A = K - \omega^2 M $$など

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

#javascript(){{
<script type="text/javascript">
  var _gaq = _gaq || [];
  _gaq.push(['_setAccount', 'UA-22313036-2']);
  _gaq.push(['_trackPageview']);
  (function() {
    var ga = document.createElement('script'); ga.type = 'text/javascript'; ga.async = true;
    ga.src = ('https:' == document.location.protocol ? 'https://ssl' : 'http://www') + '.google-analytics.com/ga.js';
    var s = document.getElementsByTagName('script')[0]; s.parentNode.insertBefore(ga, s);
  })();
</script>
}}

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

*定数行列の準備
${\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] $
//添字の連続を避けるため、ある要素qにおけるn番節点のi座標座標を指定するときは $ x_i(e_{nq}) $と呼ぶものとする。

*自由度・荷重
拘束・外力込みで定義しておく、このページでは節点における自由度を先にカウントする。
自由度に対する時間微分などあれば同様。
$ 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: は並列化推奨ループを示す。

**&bold(){parfor q = 1..q}(各要素に対するループ)
-$N^q(n, s), D^{q\prime}(d,n,s), G(d,d,s), detJ(s)$: メモリの動的確保
**&bold(){parfor s = 1..n}(数値積分各要素に関するループ)
+$\xi_{ks} = S_{ks}/\sqrt{3} $　　
+$N_n = \prod_k \frac{1 + S_{kn} \xi_{ks}}{2} $
+$D^\prime_{jn} = \frac{\partial N_n}{\partial \xi_j} =  S_{nj} \prod_k \frac{1 + \bar{\delta}_{jk} S_{nk} \xi_{ks}}{2}$
+$J_{ij} = x_{in} D^\prime_{jn} $ 　　[loop: n]
+$G_{ij} = J_{ij}^{-1} $
+$dV = det J_{ij} $
**&bold(){end }(index s)
**&bold(){end }(要素 q)

並列化パラメータの関係で一旦要素ループをここで切る。
**&bold(){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フィル
**&bold(){parfor m = 1..Σn, n = 1..Σn}(K行列等各要素に対するループ)
-$D^{n,m}(d), {\cal B}^{n,m}(d,d,d)$行列のメモリ確保
-&bold(){for s = 1..n}(数値積分): リダクション(1メモリへの複数回アクセス)が噛むので並列化しないのが安泰
+$D_{im}, D_{in} = N_{n,i} =  D^\prime_{jn}  G_{ij} $　　[loop: j]
+${\cal B}_{uvjm}, {\cal B}_{klin} = \frac{1}{2} \left[ \delta_{ik} D_{ln}  + \delta_{il} D_{kn} \right]$
+$K^q_{mn} += D_{im} D_{in} dV$　　[loop: i]
+${\cal K}^q_{ijmn} += {\cal B}_{uvjm}  {\cal C}_{uvkl} {\cal B}_{klin}  dV $ [loop: k,l,u,v ] 
+$M^q_{mn} += N_m N_n dV$
+(意図してN行列ベースの質量行列を作成したければ)${\cal M}^q_{ijmn} += \rho \delta_{ij}  N_m N_n dV$
-&bold(){end } (index s)
**&bold(){end }(index m, n)
**&bold(){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;
**&bold(){for m = 1..Lp}(各節点に対するループ)
**&bold(){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
-&bold(){for a = 1..max L(p,a)}(隣接節点に対するループ, $a > L(p,a)$ ならリターン )
-&bold(){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
-&bold(){end }(index j)
-&bold(){end }(隣接節点 a)
**&bold(){end }(index i)
**&bold(){end }(節点 m)


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

**&bold(){parfor m = 1..Lp}(各節点に対するループ)
**&bold(){parfor i = 1..Ld}(節点mの自由度)
+$r = i + L_d m$
+$l = free(r)$
+$\tilde{f}(idx_f) = f(r)   \hspace{2em} if\ l \neq -1$
-&bold(){for a = 1..max L(p,a)}(隣接節点に対するループ, $a > L(p,a)$ ならリターン )
-&bold(){for j = 1..Ld}(節点aの自由度)
+$n = m.adjacent(a)$ : 隣接節点の番号
+$qlist =m.elem(a)$ : 共有要素リスト
+$c = j + L_d a$
+$k = idx(r,c)$
+$row(k) = r   \hspace{2em} if\ k \neq -1$
+$col(k) = j + L_d n   \hspace{2em} if\ k \neq -1 $
-&bold(){for q = 1..qlist}(各要素に対するループ)
+$ K(k) += {\cal K}^q_{ijmn}   \hspace{2em} if\ k \neq -1$ 
+必要なら $M(k) += {\cal M}^q_{ijmn}  \hspace{2em} if\ k \neq -1 $
+$\tilde{f}(l) -=  {\cal K}^q_{ijmn} u( j + L_d n)  \hspace{2em} if\ l \neq -1 $
-&bold(){end }(要素 q)
-&bold(){end }(index j)
-&bold(){end }(隣接節点 a)
**&bold(){end }(index i)
**&bold(){end }(節点 m)


**最終的な連立方程式に
$ A = K - \omega^2 M $など

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

#javascript(){{
<script type="text/x-mathjax-config"> MathJax.Hub.Config({
    tex2jax: {inlineMath: [ ['$','$'], ["\\(","\\)"] ], processEscapes: true },
    TeX: { Macros: {L: '\\Large\\displaystyle'}}});  </script>
<script type="text/javascript" src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS_HTML"></script>
}}
#javascript(){{
<script type="text/javascript">
  var _gaq = _gaq || [];
  _gaq.push(['_setAccount', 'UA-22313036-2']);
  _gaq.push(['_trackPageview']);
  (function() {
    var ga = document.createElement('script'); ga.type = 'text/javascript'; ga.async = true;
    ga.src = ('https:' == document.location.protocol ? 'https://ssl' : 'http://www') + '.google-analytics.com/ga.js';
    var s = document.getElementsByTagName('script')[0]; s.parentNode.insertBefore(ga, s);
  })();
</script>
}}