三次元弾性解析コード (/3) プログラムの概要 2 年夏学期 中島研吾 科学技術計算 Ⅰ(482-27) コンピュータ科学特別講義 Ⅰ(48-24)
FEM3D-Part 2 対象とする問題 弾性体 Z ヤング率 E ポアソン比 ν U Z Z Z Y Y 直方体 一辺長さの立方体 ( 六面体 ) 要素 各方向に Y Z 個 境界条件 対称条件 U @ U Y @Y U Z @Z 強制変位 U Z @ZZ ma move
FEM3D-Part 3 構成式 : 応力 ~ひずみ関係 ヤング率 E 応力とひずみは比例 比例定数をヤング率 Eとする ( 各物質に固有の値 ) σ σ Eε ε E ポアソン比 ν 方向に荷重をかけると 横方向 (YZ) (YZ) にも変形 縮み割合をポアソン比 ν とする 各物質に固有の値 ε 金属では.3 程度 水 :.5 ゴム : ほぼ.5 非圧縮 νε σ ν E ε νε ε σ
FEM3D-Part 4 YZν.3 とすると Z U Z Z ε ε. ε νε.3 Z u u ε. 3 Y Y Y Y
FEM3D-Part 5 有限要素法の処理 支配方程式 ガラーキン法 : 弱形式 要素単位の積分 要素マトリクス生成 全体マトリクス生成 境界条件適用 連立一次方程式
FEM3D-Part 6 初期化 有限要素法の処理 : プログラム 制御変数読み込み 座標読み込み 要素生成 (: 節点数 ICELTOT: 要素数 ) 配列初期化 ( 全体マトリクス 要素マトリクス ) 要素 全体マトリクスマッピング (IndeItem) マトリクス生成 要素単位の処理 (do cel ICELTOT) 要素マトリクス計算 全体マトリクスへの重ね合わせ 境界条件の処理 連立一次方程式 共役勾配法 (CG) 応力計算
FEM3D-Part 7 三次元要素の定式化 三次元弾性力学方程式 ガラーキン法 要素マトリクス生成 宿題 プログラムの実行 データ構造 プログラムの構成 計算効率について McroAVS による可視化
FEM3D-Part 8 二次元への拡張 : 三角形要素 任意の形状を扱うことができる 特に一次要素は精度が悪く 一部の問題を除いてあまり使用されない
FEM3D-Part 9 二次元への拡張 : 四角形要素 一次元要素と同じ形状関数を 軸に適用することによって 四角形要素の定式化は可能である 三角形と比較して特に低次要素の精度はよい しかしながら 各辺が座標軸に平行な長方形でなければならない 差分法と変わらない 4 3 このような形状を扱うことができない 2
FEM3D-Part アイソパラメトリック要素 (/3) 各要素を 自然座標系 (ξη) の正方形要素 [± ±] に変換する 4 3 2 4 η 3 - - 2 各要素の全体座標系 (global coordnate)() における座標成分を 自然座標系における形状関数 [] ( 従属変数の内挿に使うのと同じ []) を使用して変換する場合 このような要素をアイソパラメトリック要素 (soparametrc element) という ξ
FEM3D-Part アイソパラメトリック要素 (2/3) 4 3 η 4 3 2 - - 2 各節点の座標 :( ) ( 2 2 ) ( 3 3 ) ( 4 4 ) ξ 各節点における 方向変位 :u u 2 u 3 u 4 u 4 ( ξ η) u 4 4 ξ η) ( ( ξ η)
FEM3D-Part 2 アイソパラメトリック要素 (3/3) η - ξ - 高次の補間関数を使えば 曲線 曲面も扱うことが可能となる そういう意味で 自然座標系 と呼んでいる
FEM3D-Part 3 2D 自然座標系の形状関数 (/2) 自然座標系における正方形上の内挿多項式は下式で与えられる : u α α 2ξ α3η α 4ξη 各節点での条件より : η 4 3 - - 2 ξ u u2 u3 u4 α α 2 4 u u2 u3 u4 α3 α 4 4 u u u2 u 4 u2 u 4 3 3 u u 4 4
FEM3D-Part 4 2D 自然座標系の形状関数 (2/2) 元の式に代入して u について整理すると以下のようになる : η 4 3 u u u u u ξ 2 2 3 3 4u4 形状関数 は以下のようになる : - 2 - ( ξ η) ( ξ )( η) 2( ξ η) ( ξ )( η) 4 4 3 ( ξ η) ( ξ )( η) 4( ξ η) ( ξ )( η) 4 4 双一次 (b-lnear) 要素とも呼ばれる 各節点における の値を計算してみよ Y 方向変位 v についても同様 v v 2v2 3v3 4v4
FEM3D-Part 5 三次元への拡張 四面体要素 : 二次元における三角形要素 任意の形状を扱うことができる 特に一次要素は精度が悪く 一部の問題を除いてあまり使用されない 高次の四面体要素は広く使用されている 本講義では低次六面体要素 ( アイソパラメトリック要素 ) を使用する tr-lnear 変位法 自由度 : 変位 各節点上で 3 成分 (uvw) が定義される
FEM3D-Part 6 3D 自然座標系の形状関数 ( ξ η ζ ) 8 2 ( ξ η ζ ) ξ η ζ 8 3 ( ξ η ζ ) ξ η ζ 8 4 ( ξ η ζ ) ξ η ζ 8 8 u ( ξ η ζ ) u 8 v ( ξ η ζ ) v 8 ( ξ )( η)( ζ ) w ( ξ η ζ ) w ( )( )( ) ( )( )( ) ( )( )( ) 5 ( ξ η ζ ) 8 6 ( ξ η ζ ) ξ η ζ 8 7 ( ξ η ζ ) ξ η ζ 8 8 ( ξ η ζ ) ξ η ζ 8 ( ) 8 7 ( ) 5 6 ( ) ( ξ η ζ ) ( ) ( ) 4 3 2 ( ξ )( η)( ζ ) ( )( )( ) ( )( )( ) ( )( )( ) ( ) ( ) ( )
FEM3D-Part 7 三次元要素の定式化 三次元弾性力学方程式 ガラーキン法 要素マトリクス生成 宿題 プログラムの実行 データ構造 プログラムの構成 計算効率について McroAVS による可視化
FEM3D-Part 8 弾性力学の支配方程式 つりあい式 ひずみ ~ 変位関係式 ひずみ ~ 応力関係式
FEM3D-Part 9 三次元のつりあい式な成応力の独立な成分は 6 つ τ τ { τ σ τ τ τ σ σ τ τ τ τ σ τ τ τ τ τ τ σ Y τ σ τ Z σ τ τ Z
FEM3D-Part 2 垂直ひずみ~ 変位の関係 PQ P Q ε u d u d ( u ) d d u d R P u / P d R v / Q Q ε ε ε u u v w
FEM3D-Part 2 せん断ひずみ~ 変位の関係 R R γ u v d u / P v / Q γ γ v w w u P d Q
FEM3D-Part 22 応力 ひずみ関係ず関係 σ σ ν ν ν ν ε ε ( ) E τ σ σ ν ν ν ν ν γ ε ε 2 ( ) ( ) ( ) E τ τ ν ν γ γ 2 2 2 ( ) τ ν γ 2
FEM3D-Part 23 ひずみ 応力関係ず関係 ν ν ν ε ε ν ν ν ν ν ν σ σ ( )( ) ( ) ( ) E γ ε ν ν ν ν τ σ 2 2 2 2 ( ) ( ) γ γ ν ν τ τ 2 2 2 2 2 [ ] D { [ ]{ ε σ D 非圧縮性材料 (ν~.5) の場合 特別な扱い必要
FEM3D-Part 24 - 方向のつりあい式に注目 σ τ τ σ τ τ ガラーキン法 T [ ] { σ τ τ dv T [ ] { [ ] { [ ] { T T T σ dv τ dv τ dv [ ] { dv
FEM3D-Part 25 - 方向のつりあい式 (/3) T T T [ ] { σ [ ] { dv σ dv [ ] { σ V Q V [ ] { ] T T σ dv [ ] { σ V S n S ds V n ds [ ] { ] T T [ ] { [ ] T σ dv σ dv { σ V 同様にして以下の弱形式が得られる ( 表面積分省略 ): V 一次元ガウスの公式 dv 部分積分の公式 T [ ] { [ ] { [ ] { T T T σ dv τ dv τ dv [ ] { dv V V V V V [ ] T { [ ] T { [ ] T T σ dv τ dv { τ dv [ ] { dv V V V
FEM3D-Part 26 V - 方向のつりあい式 (2/3) [ ] T { [ ] T { [ ] T T dv τ dv { τ dv [ ] { dv (*) σ V V V σ E ν ε νε ν 2ν νε ( )( ) ( ) [ ] E ν u νv w ν 2ν ν ( )( ) ( ) [ ] τ τ E 2 γ E [ ] u v ( ν ) 2( ν ) E E 2 γ [ ] u w ( ν ) 2( ν )
FEM3D-Part 27 要素内の変位分布形状関数要素内の変位分布 形状関数 [ ] [ ] [ ] { { { W w V v U u [ ] [ ] [ ] { { { W w V v U u [ ] [ ] [ ] { { { U u U u U u [ ] [ ] [ ] [ ] [ ] { { { { { { { V v V v U u U u U u [ ] [ ] { { W w W w E ( )( ) ( )[ ] [ ] [ ] ( ) { { { 2 E E W V U E ν ν ν ν ν σ ( ) ( ) [ ] [ ] ( ) { { 2 2 V U E E ν γ ν τ ( ) ( ) [ ] [ ] ( ) { { 2 2 W U E E ν γ ν τ
FEM3D-Part 28 - 方向のつりあい式 (3/3) 方向のつりあい式 (3/3) ( ) ν ν E E ( )E ( )( ) ( )( ) ( ) ν ν ν ν ν ν ν 2 2 2 E b E a E D とすると [ ] [ ] [ ] [ ] [ ] { { { { { V b U b W a V a U D τ σ [ ] [ ] { { W b U b τ ( ) * ( ) [ ] [ ] [ ] [ ] [ ] [ ] ( ) { { * V T T T U dv b D [ ] [ ] [ ] [ ] ( ) { { [ ] [ ] [ ] [ ] { { V T T V T T V W dv b a V dv b a [ ] { V T dv
FEM3D-Part 29 Y- 方向のつりあい式 [ ][ ] [ ] [ ] [ ] [ ] ( ) { { T T T T V T T T V dv b D [ ] [ ] [ ] [ ] ( ) { { [ ] [ ] [ ] [ ] { { [ ] { T V T T V T T dv Y W dv b a U dv b a [ ] { V dv Y Z- 方向のつりあい式 [ ][ ] [ ] [ ] [ ] [ ] ( ) { { T T T T V T T T V dv b U dv b W dv b D [ ] [ ] [ ] [ ] ( ) { { [ ] [ ] [ ] [ ] { { [ ] { T V T T V T T dv Z V dv b a U dv b a [ ] { V dv Z
FEM3D-Part 3 3 つのつりあい式 変位の 3 成分を未知数 33 つの方程式 3つの方程式はカップルしている ( 独立では無い ) 形状関数ベクトル []:8 行列 [] T [] [ ] T [ ] 等 :8 8 行列 σ τ τ τ τ σ τ τ σ Y Z
FEM3D-Part 3 要素マトリクス :24 24 行列 各節点上の (uvw) 成分が物理的にも強くカップルしているので3 自由度をまとめて扱う :8 8 行列 j ( ) 8 7 ( ) ( ) 5 6 ) ( ) 4 3 ( ) 2 ( ξ η ζ ) ( ) ( ) j j 2 j3 j2 j 22 j 32 j3 j 23 j 33 ( j K8)
FEM3D-Part 32 要素マトリクス :24 24 行列 各節点上の (uvw) 成分を独立に扱うことも可能であるが 8 8とする利点については来週以降も説明する j ( ) 8 7 ( ) ( ) 5 6 ) ( ) 4 3 ( ) 2 ( ξ η ζ ) ( ) ( ) j j 2 j3 j2 j 22 j 32 j3 j 23 j 33 ( j K8)
FEM3D-Part 33 要素マトリクス :-j j 成分 方向 (/3) j j 2 j3 j2 j 22 j 32 j3 j 23 j 33 ( j K8) ( ) 8 7 ( ) 4 3 ( ) ( ) 5 6 ( ) 2 ( ξ η ζ ) ( ) ( ) ( ) { ( ) j D j b j j V { a j b j j2 V { a j b j j 3 V dv dv dv
FEM3D-Part 34 要素マトリクス :-j j 成分 Y 方向 (2/3) j j 2 j3 j2 j 22 j 32 j3 j 23 j 33 ( j K8) ( ) 8 7 ( ) 4 3 ( ) ( ) 5 6 ( ) 2 ( ξ η ζ ) ( ) ( ) ( ) { j a j b 2 j j 22 V dv { D ( ) j b j j V { a j b j j 23 V dv dv
FEM3D-Part 35 要素マトリクス :-j j 成分 Z 方向 (3/3) j j 2 j3 j2 j 22 j 32 j3 j 23 j 33 ( j K8) ( ) 8 7 ( ) 4 3 ( ) ( ) 5 6 ( ) 2 ( ξ η ζ ) ( ) ( ) ( ) { j a j b 3 j V { a j b j j32 j33 V dv dv { D ( ) j b j j V dv
FEM3D-Part 36 あとは積分を求めれば良いあ積分を求良 自然座標系 (ξηζ) 上で定義 ガウス積分公式が ξηζ 使える ( 三次元 ) しかし 微分が M L d d d f I ) ( ζ η ξ ζ η ξ [ ] j j M j L f W W W ) ( ζ η ξ ) ( ζ ξ LM: ξηζ 方向の積分点数積分点座標値 j W W W ) ( j ζ η ξ : 積分点での重み係数 : 積分点の座標値
FEM3D-Part 37 自然座標系における偏微分 (/4) 自然座標系における偏微分 (/4) 偏微分の公式より以下のようになる : ξ ξ ξ ξ ζ η ξ ) ( η η η η ζ η ξ ξ ξ ξ ξ ) ( ζ ζ ζ ζ ζ η ξ η η η η ) ( ζ ζ ζ ζ は定義より簡単に求められるが ζ η ξ は定義より簡単に求められるが を実際の計算で使用する
FEM3D-Part 38 自然座標系における偏微分 (2/4) 自然座標系における偏微分 (2/4) マトリックス表示すると : [ ] J ξ ξ ξ ξ [ ] ζ ζ ζ η η η ζ η ζ ζ ζ ζ [ ] J : ヤコビのマトリクス [ ] J : ヤコビのマトリクス (Jacob matr Jacoban)
FEM3D-Part 39 自然座標系における偏微分 (3/4) 自然座標系における偏微分 (3/4) の定義より簡単に求められる J J 8 8 2 8 8 ξ ξ ξ ξ ξ ξ J 8 8 3 ξ ξ ξ J J 8 8 22 8 8 2 η η η η η η J 8 8 23 η η η J J 8 8 32 8 8 3 ζ ζ ζ ζ ζ ζ J 8 8 33 ζ ζ ζ
FEM3D-Part 4 自然座標系における偏微分 (4/4) 自然座標系における偏微分 (4/4) 従って下記のように偏微分を計算できる ヤコビアン (3 3 行列 ) の逆行列を求める ξ ξ ξ ξ ξ [ ] η ξ η ξ η η η ξ ξ ξ J ζ ζ ζ ζ ζ
FEM3D-Part 4 要素単位での積分 j j 2 j3 j2 j 22 j 32 j3 j 23 j 33 ( j K8) ( ) 8 7 ( ) 4 3 ( ) ( ) 5 6 ( ) 2 ( ξ η ζ ) ( ) ( ) ( ) { ( ) j D j b j j V dv V j j j D b dv
FEM3D-Part 42 自然座標系での積分 j j 2 j3 j2 j 22 j 32 j3 j 23 j 33 ( j K8) ( ) 8 7 ( ) 4 3 ( ) ( ) 5 6 ( ) 2 ( ξ η ζ ) ( ) ( ) ( ) j j j D b dv V j j j D b ddd D j j j b dt det J d ξ dηdζζ
FEM3D-Part 43 ガウスの積分公式 j j 2 j3 j2 j 22 j 32 j3 j 23 j 33 ( j K8) ( ) 8 7 ( ) 4 3 ( ) ( ) 5 6 ( ) 2 ( ξ η ζ ) ( ) ( ) ( ) D j b j j det J dξdηdζ I L M j f ( ξ η ζ ) dξdηdζ [ W ] W j W f ( ξ η j ζ )
FEM3D-Part 44 残りの手順 ここまでで 要素ごとの積分が可能となる あとは : 全体マトリクスへの重ね合わせ 境界条件処理 連立一次方程式を解く 詳細は来週以降の講義でプログラムの内容を解説 詳細は来週以降の講義でプログラムの内容を解説しながら説明する
FEM3D-Part 45 要素 全体マトリクス重ね合わせ要素全体トリク重ね合わ 5 6 4 3 (2) (2) ) ( ) ( ) ( ) ( (2) (2) (2) { ]{ [ f φ 2 5 6 2 (2) (2) 3 (2) 2 (2) (2) (2) 3 (2) 2 (2) (2) (2) (2) (2) (2) 34 (2) 33 (2) 32 (2) 3 (2) 24 (2) 23 (2) 22 (2) 2 (2) 4 (2) 3 (2) 2 (2) f f f φ φ φ 4 3 3 4 2 () () () () () () () () () { ]{ [ f f φ φ (2) 4 (2) 4 (2) 44 (2) 43 (2) 42 (2) 4 f φ 2 2 2 () () 3 () 2 () () 3 () 2 () 44 () 43 () 42 () 4 () 34 () 33 () 32 () 3 () 24 () 23 () 22 () 2 4 3 2 f f f f φ φ φ φ 2 4 4 44 43 42 4 f φ Φ Φ { ]{ [ B D F K Φ Φ Φ Φ 3 2 3 2 3 2 B B B B D D D D Φ Φ Φ 6 5 4 6 5 4 6 5 4 B B B D D D
FEM3D-Part 46 三次元要素の定式化 三次元弾性力学方程式 ガラーキン法 要素マトリクス生成 宿題 プログラムの実行 データ構造 プログラムの構成 計算効率について McroAVS による可視化
FEM3D-Part 47 宿題 ガウスの積分公式を使用して以下の四角形の面積を求めよ ( プログラムを作って 計算機で計算してください ) 4 V 3 2 :(..) 2:(4. 2.) 3:(3. 5.) 4:(2. 4.) I V dv det J dξdζξ ζ
FEM3D-Part 48 ヒント (/2) 座標値によってヤコビアン ( ヤコビの行列 ) を計算 ガウスの積分公式 (n2) に代入する I m n f j ξ η ) d ξ d η [ W ] W j f ( ξ j ) ( η mplct REAL*8 (A-HO-Z) real*8 W(2) real*8 POI(2) W().d W(2).d POI() -.577352692d POI(2).577352692d SUM.d do jp 2 do p 2 FC F(POI(p)POI(jp)) SUM SUM W (p)*w (jp)*fc enddo enddo
FEM3D-Part 49 ヒント (2/2) ξ ξ [ ] η ξ η ξ η η ξ ξ J J det η η 4 4 4 4 4 4 4 4 ξ ξ ξ ξ ξ ξ η η η η η η ( )( ) ( )( ) η ξ η ξ η ξ η ξ 4 ) ( 4 ) ( 2 ( )( ) ( )( ) η ξ η ξ η ξ η ξ 4 ) ( 4 ) ( 4 3
FEM3D-Part 5 三次元要素の定式化 三次元弾性力学方程式 ガラーキン法 要素マトリクス生成 宿題 プログラムの実行 データ構造 プログラムの構成 計算効率について McroAVS による可視化
FEM3D-Part 5 対象とする問題 弾性体 Z ヤング率 E ポアソン比 ν U Z Z Z Y Y 直方体 一辺長さの立方体 ( 六面体 ) 要素 各方向に Y Z 個 境界条件 対称条件 U @ U Y @Y U Z @Z 強制変位 U Z @ZZ ma move
FEM3D-Part 52 YZν.3 とすると Z U Z Z ε ε. ε νε.3 Z u u ε. 3 Y Y Y Y
FEM3D-Part 53 ファイルコピー インストール (/2) 三次元弾性解析コード >$ cd <$fem> >$ cp /home3/sengon/documents/class/fem/fem3d.tar. >$ tar vf fem3d.tar >$ cd fem3d >$ ls c f run FEM インストール (C) >$ cd <$fem>/fem3d/c >$ mae >$ ls../run/sol sol FEM インストール (FORTRA) >$ cd <$fem>/fem3d/f >$ mae >$ ls../run/sol sol
FEM3D-Part 54 ファイルコピー インストール (2/2) メッシュジェネレータインストール >$ cd <$fem>/fem3d/run >$ g95 O3 mgcube.f o mgcube
FEM3D-Part 55 計算の流れメッシュ生成 計算 ファイル名称固定 mgcube メッシュジェネレータ cube. メッシュファイル sol FEM 本体 IPUT.DAT 制御データ test.np 可視化用出力変位 3 成分 σ
FEM3D-Part 56 メッシュ生成 Z U Z >$ cd <$fem>/fem3d/run / >$./mgcube Y Z 各辺長さを訊いてくる このように 入れてみる Z >$ ls cube. 生成を確認 Y Y cube.
FEM3D-Part 57 制御ファイル :IPUT.DAT IPUT.DAT cube. fname METHOD PRECOD terprema( 不使用 ) 2 ITER..3 ELAST POISSO fname: メッシュファイル名 METHOD: 反復解法 ( に固定 ) PRECOD: 前処理手法 :Bloc-LU-GS:Bloc 対角スケーリング terprema: ( 不使用 ) ITER: ELAST: 反復回数上限 ヤング率 POISSO: ポアソン比 あとで.4999などの場合を試してみよ
FEM3D-Part 58 >$ cd <$fem>/fem3d/run >$./sol ( 反復法の収束履歴 ) 33 2.28867E-8 34.32592E-8 35 7.38434E-9 実行 ### DISPLACEMET at (maymazma)) 33-3.E- -3.2E-.E >$ ls test.np 生成を確認 test.np Z U Z この点での変位が表示されている 33 番節点 ( 3 ) Z Y Y
FEM3D-Part 59 三次元要素の定式化 三次元弾性力学方程式 ガラーキン法 要素マトリクス生成 宿題 プログラムの実行 データ構造 プログラムの構成 計算効率について McroAVS による可視化
FEM3D-Part 6 メッシュファイル構成 :cube. 番号は から始まっている 節点データ 節点数 節点番号 座標 要素データ 要素数 要素タイプ 要素番号 材料番号 コネクティビティ 節点グループデータ グループ数 グループ内節点数 グループ名 グループ内節点
FEM3D-Part 6 メッシュ生成コード :mgcube.f(/5) mplct REAL*8 (A-HO-Z) real(nd8) dmenson(::) allocatable :: Y real(nd8) dmenson(::) allocatable :: Y character(len8) :: GRIDFILE HHH nteger dmenson(::) allocatable :: IW!C!C -------!C IIT.!C -------!C wrte (**) 'YZ' read (**) YZ P YP Y ZP Z D.d IODTOT P*YP*ZP ICELTOT *Y *Z IBODTOT P*YP allocate (IW(IODTOT4)) IW 方向節点数 Y 方向節点数 Z 方向節点数 総節点数総要素数 Y 平面の節点数
FEM3D-Part 62 メッシュ生成コード :mgcube.f(2/5)!c cou b do ZP do j YP cou cou (-)*IBODTOT (j-)*p IW(coub) enddo enddo cou b 2 do ZP j do P cou cou (-)*IBODTOT (j-)*p IW(coub) enddo enddo cou b 3 do j YP do P cou cou (-)*IBODTOT (j-)*p IW(coub) enddo enddo cou b 4 ZP do j YP do P cou cou (-)*IBODTOT (j-)*p IW(coub) enddo enddo mn の節点を IW(b) に格納 (byp*zp) Z Y U Z Z Y
FEM3D-Part 63 メッシュ生成コード :mgcube.f(2/5)!c cou b do ZP do j YP cou cou (-)*IBODTOT (j-)*p IW(coub) enddo enddo cou b 2 do ZP j do P cou cou (-)*IBODTOT (j-)*p IW(coub) enddo enddo cou b 3 do j YP do P cou cou (-)*IBODTOT (j-)*p IW(coub) enddo enddo cou b 4 ZP do j YP do P cou cou (-)*IBODTOT (j-)*p IW(coub) enddo enddo YYmnの節点をIW(b2) に格納 (bp*zp) ZP) Z U Z Z Y Y
FEM3D-Part 64 メッシュ生成コード :mgcube.f(2/5)!c cou b do ZP do j YP cou cou (-)*IBODTOT (j-)*p IW(coub) enddo enddo Z U Z cou b 2 do ZP Z j do P cou cou (-)*IBODTOT (j-)*p IW(coub) Y enddo enddo cou Y b 3 do j YP do P cou cou (-)*IBODTOT (j-)*p IW(coub) enddo enddo cou b 4 ZP do j YP do P cou cou (-)*IBODTOT (j-)*p IW(coub) enddo enddo ZZmn の節点を IW(b3) に格納 (bp*yp) ZZmaの節点をIW(b4) に格納 (bp*yp) YP)
FEM3D-Part 65 メッシュ生成コード :mgcube.f(3/5)!c!c -------------!C GeoFEM data!c -------------!C wrte (**) 'GeoFEM grdfle name?' GRIDFILE 'cube.' open (2 fle GRIDFILE status'unnown'form'formatted') wrte(2'()') IODTOT cou do ZP do j YP do P dfloat(-)*d YY dfloat(j-)*d ZZ dfloat(-)*d cou cou wrte (2'(3(pe6.6))') cou YY ZZ enddo enddo enddo wrte(2'()') ICELTOT IELMTYPL 36 wrte(2'()') (IELMTYPL ICELTOT) 節点数節点番号 座標 要素数要素タイプ ( 使わないデータ )
FEM3D-Part 66 cube.: 節点データ 要素数 要素タイプ (YZ4) 25 5*5*5.E.E.E 2.E.E.E 3 2.E.E.E 4 3.E.E.E 5 4.E.E.E 6.E.E.E 7.E.E.E 8 2.E.E.E 9 3.E.E.E ( 途中省略 ) 2.EE 4.E 4.E 22.E 4.E 4.E 23 2.E 4.E 4.E 24 3.E 4.E 4.E 25 4.E 4.E 4.E 64 4*4*4 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 36 move
FEM3D-Part 67 メッシュ生成コード :mgcube.f(4/5) cou mat do Z do j Y do cou cou n (-)*IBODTOT (j-)*p n2 n n3 n2 P n4 n3 - n5 n IBODTOT n6 n2 IBODTOT n7 n3 IBODTOT n8 n4 IBODTOT wrte (2'()') cou mat n n2 n3 n4 & & n5 n6 n7 n8 enddo enddo enddo mat: 材料番号 () 8 7 5 6 4 3 2
FEM3D-Part 68 cube.: 要素データ 2 7 6 26 27 32 3 2 2 3 8 7 27 28 33 32 3 3 4 9 8 28 29 34 33 4 4 5 9 29 3 35 34 5 6 7 2 3 32 37 36 6 7 8 3 2 32 33 38 37 7 8 9 4 3 33 34 39 38 8 9 5 4 34 35 4 39 9 2 7 6 36 37 42 4 2 3 8 7 37 38 43 42 3 4 9 8 38 39 44 43 2 4 5 2 9 39 4 45 44 3 6 7 22 2 4 42 47 46 ( 途中省略 ) 42 62 63 68 67 87 88 93 92 43 63 64 69 68 88 89 94 93 44 64 65 7 69 89 9 95 94 45 66 67 72 7 9 92 97 96 46 67 68 73 72 92 93 98 97 47 68 69 74 73 93 94 99 98 48 69 7 75 74 94 95 99 49 76 77 82 8 2 7 6 5 77 78 83 82 2 3 8 7 5 78 79 84 83 3 4 9 8 52 79 8 85 84 4 5 9 53 8 82 87 86 6 7 2 54 82 83 88 87 7 8 3 2 55 83 84 89 88 8 9 4 3 56 84 85 9 89 9 5 4 57 86 87 92 9 2 7 6 58 87 88 93 92 2 3 8 7 59 88 89 94 93 3 4 9 8 6 89 9 95 94 4 5 2 9 6 9 92 97 96 6 7 22 2 62 92 93 98 97 7 8 23 22 63 93 94 99 98 8 9 24 23 64 94 95 99 9 2 25 24
FEM3D-Part 69 メッシュ生成コード :mgcube.f(5/5) IGTOT 4 IBT YP*ZP IBT2 P*ZP IBT IBT3 P*YP IBT2 IBT4 P*YP IBT3 wrte (2'()') IGTOT wrte (2'()') IBT IBT2 IBT3 IBT4 HHH 'mn' wrte (2'(a8)') HHH wrte (2'()') (IW() YP*ZP) HHH 'Ymn' wrte (2'(a8)') HHH wrte (2'()') ( ) (IW(2) ( P*ZP) HHH 'Zmn' wrte (2'(a8)') HHH wrte (2'()') (IW(3) P*YP) HHH 'Zma' wrte (2'(a8)') HHH wrte (2'()') (IW(4) P*YP) IGTOT IBT グループ総数 (mnymnzmnzma) Zma) 累積数 ( 以下略 )!C stop end deallocate (IW) close (2)
FEM3D-Part 7 cube.: 節点グループデータ mn Ymn Zmn Zma 4 25 5 75 6 6 2 26 3 36 4 46 5 56 6 66 7 76 8 86 9 96 6 6 2 2 3 4 5 26 27 28 29 3 5 52 53 54 55 76 77 78 79 8 2 3 4 5 2 3 4 5 6 7 8 9 2 3 4 5 6 7 8 9 2 2 22 23 24 25 2 3 4 5 6 7 8 9 2 3 4 5 6 7 8 9 2 2 22 23 24 25 ( 以下使用せず )
FEM3D-Part 7 メッシュ生成 実は技術的には大きな課題 複雑形状 大規模メッシュ 並列化が難しい 市販のメッシュ生成アプリケーション FEMAP CAD データとのインタフェース move
FEM3D-Part 72 三次元要素の定式化 三次元弾性力学方程式 ガラーキン法 要素マトリクス生成 宿題 プログラムの実行 データ構造 プログラムの構成 計算効率について McroAVS による可視化
FEM3D-Part 73 有限要素法の処理 : プログラム 初期化 制御変数読み込み 座標読み込み 要素生成 (: 節点数 E: 要素数 ) 配列初期化 ( 全体マトリクス 要素マトリクス ) 要素 全体マトリクスマッピング (IndeItem) マトリクス生成 要素単位の処理 (do cel E) 要素マトリクス計算 全体マトリクスへの重ね合わせ 境界条件の処理 連立一次方程式 共役勾配法 (CG) 応力計算
FEM3D-Part 74 test メインプログラム nput_cntl 制御データ入力 nput_grd メッシュファイル入力 fnd_node 節点探索 mat_con 行列コネクティビティ生成 msort ソート mat_con 行列コネクティビティ生成 mat_ass_manass man 係数行列生成 jacob ヤコビアン計算 mat_ass_bc 境界条件処理 三次元弾性解析コードfem3d の構成 solve33 線形ソルバー制御 recover_stress 応力計算 output_ucd 可視化処理 cg_3 CG 法計算 jacob ヤコビアン計算
FEM3D-Part 75 全体処理 #nclude <stdo.h> #nclude <stdlb.h> FILE* fp_log; #defne GLOBAL_VALUE_DEFIE #nclude "pfem_utl.h" etern vod IPUT_CTL(); etern vod IPUT_GRID(); etern vod MAT_CO(); etern vod MAT_CO(); etern vod MAT_ASS_MAI(); etern vod MAT_ASS_BC(); ASS etern vod SOLVE33(); etern vod RECOVER_STRESS(); etern vod OUTPUT_UCD(); nt man() { /** Logfle for debug **/ f( (fp_logfopen("log.log""w")) ULL){ fprntf(stdout"nput fle cannot be opened! n"); et(); IPUT_CTL(); IPUT_GRID(); MAT_CO(); MAT_CO(); MAT_ASS_MAI(); MAT_ASS_BC() ; SOLVE33(); RECOVER_STRESS(); OUTPUT_UCD() ;
FEM3D-Part 76 Global 変数表 :pfem_utl.h/f(/3) 変数名種別サイズ I/O 内容 fname C [8] I メッシュファイル名 P I I 節点数 ICELTOT I I 要素数 ODGRPtot I I 節点グループ数 YZ R [][3] I 節点座標 ICELOD I [ICELTOT][8] I 要素コネクティビティ ODGRP_IDE I [ODGRPtot] I 各節点グループに含まれる節点数 ( 累積 ) ODGRP_ITEM I [ODGRP_IDE[ODG RPTOT]] I 節点グループに含まれる節点 ODGRP_AME C8 [ODGRP_IDE[ODG RPTOT]] I 節点グループ名 L U I O 各節点非対角成分数 ( 上三角 下三角 ) PL PU I O 非対角成分総数 ( 上三角 下三角 ) D R [9*] O 全体行列 : 対角ブロック B R [3*] O 右辺ベクトル 未知数ベクトル
FEM3D-Part 77 Global 変数表 :pfem_utl.h/f(2/3) 変数名種別サイズ I/O 内容 ALUG R [9*] O 全体行列 : 対角ブロックの完全 LU 分解 ALAU R [9*PL][9*PU] O 全体行列 : 上 下三角成分 ndel ndeu I [] O 全体行列 : 非零非対角ブロック数 teml temu I [PL] [PU] O 全体行列 : 上 下三角ブロック ( 列番号 ) IL IU I [] O 各節点の上 下三角ブロック数 IALIAUIAU I [][L][][U][][U] O 各節点の上 下三角ブロック ( 列番号 ) IWK I [][2] O ワーク用配列 METHOD I I 反復解法 ( に固定 ) PRECOD I I 前処理手法 (: ブロック SSOR: ブロック対角スケーリング ) ITER ITERactual I I 反復回数の上限 実際の反復回数 RESID R I 打ち切り誤差 (.e-8 に設定 ) SIGMA_DIAG R I LU 分解時の対角成分係数 (. に設定 ) pfemiarra I [] O 諸定数 ( 整数 ) pfemrarra R [] O 諸定数 ( 実数 )
FEM3D-Part 78 Global 変数表 :pfem_utl.h/f(3/3) 変数名種別サイズ I/O 内容 O8th R I.25 PQ PE PT R [2][2][8] O 各ガウス積分点における ξ η ζ ( ~ 8) POS WEI R [2] O 各ガウス積分点の座標 重み係数 COL COL2 I [] O ソート用ワーク配列 SHAPE R [2][2][2][8] O 各ガウス積分点における形状関数 (~8) P PY PZ R [2][2][2][8] O 各ガウス積分点における ( ~ 8 ) DETJ R [2][2][2] O 各ガウス積分点におけるヤコビアン行列式 ELAST POISSO R I ヤング率 ポアソン比 SIGMA_ TAU_ R [][3] O 節点における垂直 せん断応力成分
FEM3D-Part 79 制御ファイル入力 :IPUT_CTL #nclude <stdo.h> #nclude <stdlb.h> #nclude "pfem_utl.h" /** **/ vod IPUT_CTL() { FILE *fp; f( (fpfopen("iput.dat""r")) ULL){ fprntf(stdout"nput fle cannot be opened! n"); et(); fscanf(fp"%s"fname); f(f " f fscanf(fp"%d %d"&method&precod); fscanf(fp"%d"&terprema); fscanf(fp"%d"&iter); fscanf(fp "%lf %lf" &ELAST &POISSO); fclose(fp); f( ( terprema < ) ){ terprema ; f( ( terprema > 4 ) ){ terprema 4; SIGMA_DIAG.; SIGMA.; RESID.e-8; SET ; pfemrarra[] RESID; pfemrarra[] SIGMA_DIAG; pfemrarra[2] SIGMA; IPUT.DAT cube. fname METHOD PRECOD 2 terprema( 不使用 ) ITER..3 ELAST POISSO pfemiarra[] ITER; pfemiarra[] METHOD; pfemiarra[2] PRECOD; pfemiarra[3] SET; pfemiarra[4] terprema;
FEM3D-Part 8 メッシュ入力 :IPUT_GRID(/2) #nclude <stdo.h> #nclude <stdlb.h> #nclude "pfem_utl.h" #nclude "allocate.h" vod IPUT_GRID() { FILE *fp; nt jnncelse; nt TYPEIMAT; /** **/ f( (fpfopen(fname"r")) ULL){ fprntf(stdout"nput f( fle cannot be opened! n"); et(); ODE fscanf(fp"%d"&); P; YZ(KREAL**)allocate_matr(seof(KREAL)3); for(;<;){ for(j;j<3;j){ YZ[][j].; for(;<;){ for( ;<;){ fscanf(fp"%d %lf %lf %lf"&&yz[][]&yz[][]&yz[][2]);
FEM3D-Part 8 allocate deallocate 関数 #nclude <stdo.h> #nclude <stdlb.h> vod* allocate_vector(nt sent m) { vod d*; *a; f ( ( a(vod * )malloc( m * se ) ) ULL ) { fprntf(stdout"error:memor does not enough! n vector n"); et(); return a; allocate を FORTRA 並みに簡単にやるための関数 vod deallocate_vector(vod *a) { free( a ); vod** allocate_matr(nt sent mnt n) { vod **aa; nt ; f((aa(vod ( (vod ** )malloc( m * seof(vod*) )) ) ULL ) { fprntf(stdout"error:memor does not enough! aa n matr n"); et(); f ( ( aa[](vod * )malloc( m * n * se ) ) ULL ) { fprntf(stdout"error:memor does not enough! n matr n"); et(); for(;<m;) aa[](char*)aa[-]se*n; return aa; vod deallocate matr(vod **aa) vod deallocate_matr(vod **aa) { free( aa );
FEM3D-Part 82 メッシュ入力 :IPUT_GRID(2/2) /** **/ ELEMET fscanf(fp"%d"&iceltot); ICELOD(KIT**)allocate_matr(seof(KIT)ICELTOT8); t( ICELTOT for(;<iceltot;) fscanf(fp"%d"&type); ICELOD[][j] の中身としては から始まる通し節点番号がそのまま読み込まれている 要素番号は から番号付け /** **/ for(cel;cel<iceltot;cel){ fscanf(fp"%d %d %d %d %d %d %d %d %d %d"&&imat &ICELOD[cel][]&ICELOD[cel][]&ICELOD[cel][2]&ICELOD[cel][3] l][] l][2] l][3] &ICELOD[cel][4]&ICELOD[cel][5]&ICELOD[cel][6]&ICELOD[cel][7]); ODE grp. nfo. fscanf(fp"%d"&odgrptot); ODGRP_IDE(KIT* )allocate_vector(seof(kit)odgrptot); ODGRP_AME (CHAR8*)allocate_vector(seof(CHAR8)ODGRPtot); for(;<odgrptot;) ODGRP_IDE[]; for(;<odgrptot;) fscanf(fp"%d"&odgrp_ide[]); nnodgrp_ide[odgrptot]; ODGRP_ITEM(KIT*)allocate_vector(seof(KIT)nn); for(;<odgrptot;){ S ODGRP_IDE[]; E ODGRP_IDE[]; fscanf(fp"%s"odgrp_ame[].name); nn E - S; f( nn! ){ for(s;<e;) fscanf(fp"%d"&odgrp_item[]); fclose(fp);
FEM3D-Part 83 メッシュ入力 :IPUT_GRID(2/2) /** **/ ELEMET fscanf(fp"%d"&iceltot); ICELOD(KIT**)allocate_matr(seof(KIT)ICELTOT8); t( ICELTOT for(;<iceltot;) fscanf(fp"%d"&type); /** **/ for(cel;cel<iceltot;cel){ fscanf(fp"%d %d %d %d %d %d %d %d %d %d"&&imat &ICELOD[cel][]&ICELOD[cel][]&ICELOD[cel][2]&ICELOD[cel][3] l][] l][2] l][3] &ICELOD[cel][4]&ICELOD[cel][5]&ICELOD[cel][6]&ICELOD[cel][7]); ODE grp. nfo. fscanf(fp"%d"&odgrptot); ODGRP_IDE(KIT* )allocate_vector(seof(kit)odgrptot); ODGRP_AME (CHAR8*)allocate_vector(seof(CHAR8)ODGRPtot); for(;<odgrptot;) ODGRP_IDE[]; for(;<odgrptot;) fscanf(fp"%d"&odgrp_ide[]); nnodgrp_ide[odgrptot]; ODGRP_ITEM(KIT*)allocate_vector(seof(KIT)nn); for(;<odgrptot;){ S ODGRP_IDE[]; E ODGRP_IDE[]; fscanf(fp"%s"odgrp_ame[].name); nn E - S; f( nn! ){ for(s;<e;) fscanf(fp"%d"&odgrp_item[]); 節点グループの中身も から始まる通し節点番号がそのまま読み込まれている fclose(fp);
FEM3D-Part 84 fem3d: いくつかの特徴 非対角成分 上三角 下三角を分けて記憶 U ブロックとして記憶 ベクトル : 節点 3 成分 L 行列 : 各ブロック 9 成分 行列の各成分ではなく 節点上の3 変数に基づくブロックとして処理する
有限要素法で得られるマトリクス FEM3D-Part 85 疎行列 Φ Φ Φ 3 2 3 2 F F F D D D 疎行列 が多い A( j) のように正方行列の Φ Φ Φ Φ 7 6 5 4 3 7 6 5 4 3 F F F F D D D D A(j) のように正方行列の全成分を記憶することは疎行列では非効率的 Φ Φ Φ Φ Φ 2 9 8 2 9 8 F F F F F D D D D D 行列では非効率的 密 行列向け Φ Φ Φ Φ Φ 6 5 4 3 2 6 5 4 3 2 F F F F F D D D D D 有限要素法 : 非零非対角成分の数は高々 数百 規模 例えば未知数が 8 個あるとすると記憶容量 ( ワード数 ) は例えば未知数が 個あるとすると記憶容量 ( ワド数 ) は 正方行列 :O( 6 ) 非零非対角成分数 :O( ) 非零成分のみ記憶するのが効率的
FEM3D-Part 86 d.fd.c におけるマトリクス関連変数 変数名型サイズ内容 I - 未知数総数 PLU I - 連立一次方程式係数マトリクス非対角成分総数 Dag(:) R 連立一次方程式係数マトリクス対角成分 U(:) R 連立一次方程式未知数ベクトル Rhs(:) R 連立一次方程式右辺ベクトル Inde(:) I : Item(:) I PLU AMat(:) R PLU 係数マトリクス非対角成分要素番号用一次元圧縮配列 ( 非対角成分数 ) 係数マトリクス非対角成分要素番号用一次元圧縮配列 ( 非対角成分要素 ( 列 ) 番号 ) 係数マトリクス非対角成分要素番号用一次元圧縮配列 ( 非対角成分 ) 非零非対角成分のみを格納する Compressed Row Storage 法を使用している
行列ベクトル積への適用 ( 非零 ) 非対角成分のみを格納疎行列向け方法 FEM3D-Part 87 ( 非零 ) 非対角成分のみを格納 疎行列向け方法 Compressed Row Storage (CRS) 対角成分 ( 実数 ) Dag () 対角成分 ( 実数 ) Inde() 非対角成分数に関する一次元配列 ( 通し番号 ) ( 整数 ) () 非対角成分の要素 ( 列 ) 番号 Item() 非対角成分の要素 ( 列 ) 番号 ( 整数 nde()) AMat() 非対角成分 ( 実数 d ()) ( 実数 nde()) {Y [A]{ Φ Φ 2 2 F F D D do Y() Dag()*() d I d ( ) I d () Φ Φ Φ Φ Φ Φ 8 7 6 5 4 3 8 7 6 5 4 3 F F F F F F D D D D D D do Inde(-) Inde() Y() Y() Amat()*(Item()) enddo Φ Φ Φ Φ Φ Φ 4 3 2 9 4 3 2 9 F F F F F F D D D D D D enddo Φ Φ 6 5 6 5 F F D D
FEM3D-Part 88 行列ベクトル積への適用 ( 非零 ) 非対角成分のみを格納 疎行列向け方法 Compressed Row Storage (CRS) {Q[A]{P for(;<;){ W[Q][] Dag[] * W[P][]; for(inde[];<inde[];){ W[Q][] AMat[]*W[P][Item[]]; []]
FEM3D-Part 89 行列ベクトル積 : 密行列 とても簡単 a a2... a a a a a a 2 22 2 2......... 2 2 M M a a a a 2 a a2... a a {Y [A]{ do j Y(j).d do Y(j) Y(j) A(j)*() enddo enddo
FEM3D-Part 9 Compressed Row Storage (CRS) 2 3 4 5 6 7 8. 2.4 3.2 2 4.3 3.6 2.5 3.7 9. 3 5.7.5 3. 4 4. 9.8 2.5 2.7 5 3. 9.5.4.5 4.3 6 6.5 2.4 9.5 7 6.4 2.5.4 23. 3. 8 9.5.3 9.6 3. 5. 3
FEM3D-Part 9 Compressed Row Storage (CRS):C 2 3 4 5 6 7 2 3 4 5 6 7. 24 2.4 32 3.2 8 4 4.3 3.6 2.5 3.7 9. 対角成分 3 5 7 Dag[]. 5.7.5 3. Dag[] 3.6 2 4 6 Dag[2] 5.7 4 4. 98 9.8 25 2.5 27 2.7 Dag[3] 9.8 3 4 5 Dag[4].5 Dag[5] 2.4 3. 9.5.4.5 4.3 Dag[6] 23. 2 4 6 Dag[7] 5.3 6.5 2.4 9.5 2 5 6 64 6.4 25 2.5 4.4 23. 3. 2 5 6 7 9.5.3 9.6 3. 5.3 2 3 5 7
FEM3D-Part 92 Compressed Row Storage (CRS):C 2 3 4 5 6 7 2 3 4 5 6 7. 24 2.4 32 3.2 4 3.6 4.3 2.5 3.7 9. 3 5 7 5.7.5 3. 2 4 6 98 9.8 4 4. 25 2.5 27 2.7 3 4 5.5 3. 9.5.4 4.3 4 2 6 2.4 6.5 9.5 5 2 6 23. 64 6.4 25 2.5 4.4 3. 6 2 5 7 5.3 9.5.3 9.6 3. 7 2 3 5
FEM3D-Part 93 Compressed Row Storage (CRS):C 非対角成分数 Inde[]. 3.6 24 2.4 4.3 32 3.2 2.5 3.7 9. 4 3 5 7 2 Inde[] 2 4 Inde[2] 6 2 5.7.5 3. 2 4 6 2 Inde[3] 8 3 98 9.8 4 4. 25 2.5 27 2.7 3 4 5 3 Inde[4] 4.5 3. 9.5.4 4.3 4 2 6 4 Inde[5] 5 5 2.4 6.5 9.5 5 2 6 2 Inde[6] 7 6 23. 64 6.4 25 2.5 4.4 3. 6 2 5 7 4 Inde[7] 2 7 5.3 9.5.3 9.6 3. 4 Inde[8] 25 PLU 25 7 2 3 5 (Inde[]) Inde[]~Inde[]- 番目が 行目の非対角成分
FEM3D-Part 94 Compressed Row Storage (CRS):C 2 3 4 5 6 7 非対角成分数 Inde[]. 24 2.4 32 3.2 4 2 Inde[] 2 3.6 4.3 2.5 3.7 9. 3 5 7 4 Inde[2] 6 5.7.5 3. 2 4 6 2 Inde[3] 8 98 9.8 4 4. 25 2.5 27 2.7 3 4 5 3 Inde[4].5 3. 9.5.4 4.3 4 2 6 4 Inde[5] 5 2.4 6.5 9.5 5 2 6 2 Inde[6] 7 23. 64 6.4 25 2.5 4.4 3. 6 2 5 7 4 Inde[7] 2 5.3 9.5.3 9.6 3. 7 2 3 5 4 Inde[8] 25 Inde[]~Inde[]- 番目が 行目の非対角成分 PLU 25 (Inde[])
FEM3D-Part 95 Compressed Row Storage (CRS):C 2 3 4 5 6 7. 24 2.4 32 3.2 4 3.6 4.3 2.5 3.7 9. 3 5 7 5.7.5 3. 2 4 6 98 9.8 4 4. 25 2.5 27 2.7 3 4 5.5 3. 9.5.4 4.3 4 2 6 2.4 6.5 9.5 5 2 6 23. 64 6.4 25 2.5 4.4 3. 6 2 5 7 5.3 9.5.3 9.6 3. 7 2 3 5 例 : Item[ 6] 4 AMat[ 6].5 Item[8] 2 AMat[8] 2.5
FEM3D-Part 96 Compressed Row Storage (CRS):C 2 3 4 5 6 7. 24 2.4 32 3.2 3.6 4.3 4 2.5 3.7 9. Dag [] 対角成分 ( 実数 []) Inde[] 非対角成分数に関する一次元配列 22 333 544 755 ( 通し番号 )( 整数 []) 5.7.5 3. Item[] 非対角成分の要素 ( 列 ) 番号 2 46 67 ( 整数 [Inde[]]) 98 9.8 4 4. 25 2.5 27 2.7 AMat[] 非対角成分 3 8 49 5 ( 実数 [Inde[]]).5 3. 9.5.4 4.3 4 22364 3 4 {Y[A]{ 2.4 6.5 9.5 for(;<;){ 5 25 66 Y[] Dag[] * []; 23. 64 6.4 25 2.5 4.4 3. for(inde[];<inde[];){ 6 7 28 5972 Y[] AMat[]*[Item[]]; 5.3 9.5.3 9.6 3. 7 2 222 323524 24
FEM3D-Part 97 ブロックとして記憶 (/3) 記憶容量が減る Inde Item に関する記憶容量を削減できる j j j j
FEM3D-Part 98 ブロックとして記憶 (2/3) 計算効率 間接参照 ( メモリに負担 ) と計算の比が大きくなる ベクトル スカラー共に効く :2 倍以上の性能 連続領域 キャッシュに載る ループあたりの計算量増加 do 3* do Y() D()*() (3*-2) do nde(-) nde() 2 (3*-) tem() 3 (3*) enddo enddo Y() Y() AMAT()*() Y(3*-2) D(9*-8)*D(9*-7)*2D(9*-6)*3 Y(3*-) D(9*-5)*D(9*-4)*2D(9*-3)*3 Y(3*I ) D(9*-2)*D(9*-)*2D(9*I )*3 do nde(-) nde() enddo enddo tem() (3*-2) 2 (3*-) 3 (3*) Y(3*-2) Y(3*-2)AMAT(9*-8)*AMAT(9*-7)*2 & AMAT(9*-6)*3 Y(3*-) Y(3*-)AMAT(9*-5)*AMAT(9*-4)*2 & AMAT(9*-3)*3 Y(3*I ) Y(3*I )AMAT(9*-2)*AMAT(9*-)*2 & AMAT(9* )*3
FEM3D-Part 99 三次元要素の定式化 三次元弾性力学方程式 ガラーキン法 要素マトリクス生成 宿題 プログラムの実行 データ構造 プログラムの構成 計算効率について McroAVS による可視化
FEM3D-Part マイクロプロセッサの動向 CPU 性能 メモリバンド幅のギャップ http://www.streambench.org/
FEM3D-Part スカラープロセッサ CPU- キャッシュ - メモリの階層構造 FAST CPU Regster Cache 小容量 (MB): 一時置き場高価大きい ( 億以上のトランジスタ ) SLOW Man Memor 大容量 (GB) 廉価
FEM3D-Part 2 ベクトルプロセッサ ベクトルレジスタと高速メモリ Vector Processor 単純構造のDOループの並列処理 単純 大規模な演算に適している Vector Regster do A() B() C() enddo Ver FAST Man Memor
FEM3D-Part 3 典型的な挙動 : CG 法 ( 疎行列 FEM) 3..E2 2.5 4 % of pea GFL LOPS 2..5..5 8 % of pea GFL LOPS.E.E..E4.E5.E6.E7 IBM-SP3: DOF: Problem Se 問題サイズが小さい場合はキャッシュの影響のため性能が良い.E-.E4.E5.E6.E7 DOF: Problem Se Earth Smulator: 大規模な問題ほどベクトル長が長くなり 性能が高い
FEM3D-Part 4 プロセッサに応じたチューニング メモリ参照の最適化 に尽きる j A(j)
FEM3D-Part 5 プロセッサに応じたチューニング ( 続き ) ベクトルプロセッサ ループ長を大きくとる スカラープロセッサ キャッシュを有効利用 細切れにしたデータを扱う PC クラスタなどでは キャッシュサイズを大きくできない一方で メモリレイテンシ ( 立ち上がりオーバーヘッド ) の減少 メモリバンド幅の増加の傾向 しかしマルチコア化で帳消し 共通事項 メモリアクセスの連続性 局所性 計算順序の変更によって計算結果が変わる可能性について注意すること
FEM3D-Part 6 Scalar スカラープロセッサの代表的なチューニング ループアンローリング ループのオーバーヘッドの削減 ロード ストアの削減 ブロック化 キャッシュミスの削減
FEM3D-Part 7 Scalar ループアンローリング ロード ストアの削減 (/4) ループ処理に対する演算の割合が増す do j do A() A() B()*C(j) enddo enddo do j - 2 do A() A() B()*C(j) A() A() B()*C(j) enddo enddo T2K 東大での計算時間 4.23438E- 3.85938E- 2.6788E- do j -3 4 do A() A() B()*C(j) A() A() B()*C(j) A() A() B()*C(j2) A() A() B()*C(j3) enddo enddo
FEM3D-Part 8 Scalar ループアンローリング ロード ストアの削減 (2/4) ロード : メモリ キャッシュ レジスタ ストア : ロードの逆 ロード ストアが少ないほど効率良い FAST CPU Regster Cache SLOW Man Memor
FEM3D-Part 9 Scalar ループアンローリング ロード ストアの削減 (3/4) do j do A() A() B()*C(j) ストア ロード ロードロード enddo enddo A()B()C(j) に対して 各ループでロード ストアが () ()C( j) に対して 各ルプでドストアが発生する
FEM3D-Part Scalar ループアンローリング ロード ストアの削減 (4/4) do j -3 4 do A() A() B()*C(j) ロードロードロード A() A() B()*C(j) A() A() B()*C(j2) A() A() B()*C(j3) ストア enddo enddo 同一ループ内で複数回表れている場合には 最初にロード 最後にストアが実施され その間レジスターに保持される 計算順序に注意
FEM3D-Part Scalar ループ入替によるメモリ参照最適化 (/2) TYPE-A do do j A(j) A(j) B(j) enddo enddo TYPE-B do j do A(j) A(j) B(j) enddo enddo j A(j) FORTRA では A(j) のアドレスは A() ) A(2) A(3) A() A(2) A(22) A() A(2) A() のように並んでいる C は逆 :A[][] A[][] A[][2] A[-][] A[-][] A[-][-] この順番にアクセスしないと効率悪い ( ベクトル スカラーに共通 )
FEM3D-Part 2 Scalar ループ入替によるメモリ参照最適化 (2/2) TYPE-A do do j A(j) A(j) B(j) enddo enddo TYPE-B do j do enddo enddo A(j) A(j) B(j) TYPE-A for (j; j<; j){ for (; <; ){ A[][j] A[][j] B[][j]; TYPE-B for (; <; ){ for (j; j<; j){ A[][j] A[][j] B[][j]; T2K 東大での計算時間 (FORTRA) ### ### 52 A 3.25E-2 B 3.9625E-3 ### ### 24 A 2.34375E- B 7.825E-3 ### ### 536 A 3.476563E- B.5625E-2 2 ### ### 248 A 9.296875E- B 3.25E-2 ### ### 256 A 9.6875E- B 4.6875E-2 ### ### 372 A 2.52344E B 7.42875E-2 ### ### 3584 A.92875E B 9.765625E-2 ### ### 496 A 3.84688E B.25E-
FEM3D-Part 3 Scalar ブロック化によるキャッシュミス削減 (/7) do do j A(j) A(j) B(j) enddo enddo for (j;j<;j){ for (;<;){ A[j][] A[j][] B[][j]; 計算の処理の都合でこのような順番で計算せざるを得ない場合
FEM3D-Part 4 Scalar キャッシュの有効利用 キャッシュ 実際は64bte~28bteの細かいキャッシュラインに分かれている キャッシュライン単位でメモリへのリクエストが行われる FAST SLOW CPU Regster Cache Man Memor 小容量 (MB) 高価大きい ( 億以上のトランジスタ ) 大容量 (GB) 廉価 TLB(Translate Looasde Buffer) アドレス変換バッファ 仮想アドレスから実アドレスへの変換機能 TLB 用のキャッシュ 通常 28 8 bte 程度 : リンク時に可変 キャッシュ が有効に使われていれば TLB ミスも起こりにくい
FEM3D-Part 5 Scalar ブロック化によるキャッシュミス削減 (2/7) AB でメモリアクセスパターンが相反 特に B は効率が悪い j j A[j][] B[][j]
FEM3D-Part 6 Scalar ブロック化によるキャッシュミス削減 (3/7) 例えば キャッシュラインサイズを 4 ワードとすると配列の値は以下のようにキャッシュに転送される j A[j][]
FEM3D-Part 7 Scalar ブロック化によるキャッシュミス削減 (4/7) したがって A[][][ ] をアクセスしたら A[][][ ] A[][] ] A[][2] ] A[][3] が A[][9] をアクセスしたらA[][9] A[2][] A[2][] A[2][2] がそれぞれキャッシュ上にあるということになる 2 j 2 2 22 3 A[j][] 9
FEM3D-Part 8 Scalar ブロック化によるキャッシュミス削減 (5/7) したがって 以下のようなブロック型のパターンでアクセスすれば キャッシュを有効利用できる j j A[j][] B[][j]
FEM3D-Part 9 Scalar ブロック化によるキャッシュミス削減 (6/7) で囲んだ部分はキャッシュに載っている j j A[j][] B[][j]
FEM3D-Part 2 Scalar ブロック化によるキャッシュミス削減 (7/7) 2 2 ブロック for (j;j<;j){ for (;<;){ A[j][] A[j][] B[][j]; for (j;j<-;j2){ for (;<-;2){ A[j ][ ] A[j ][ ] B[ ][j ]; A[j ][] A[j ][] B[ ][j]; A[j][ ] A[j][ ] B[][j ]; A[j][] A[j][] B[][j]; T2K 東大での実行時間 (FORTRA) ### ### 24 BASIC 2.734375E-2 22.5625E-2 ### ### 536 BASIC 6.25E-2 22 3.55625E-2 ### ### 372 BASIC 2.57825E- 22.484375E- ### ### 3584 BASIC 3.7938E- 22 2.325E- ### ### 496 BASIC 8.4375E- 22 4.375E-
FEM3D-Part 2 Scalar チューニング : まとめ スカラープロセッサ 密行列 疎行列の場合はもっと難しい ( 研究途上の課題 ) しかし 基本的な考え方は変わらない メモリアクセスの効率化
FEM3D-Part 22 ブロックとして記憶 (3/3) 計算の安定化 対角成分で割るのではなく 対角ブロックの完全 LU 分解を求めて解く 特に悪条件問題で有効 j j j j
FEM3D-Part 23 三次元要素の定式化 三次元弾性力学方程式 ガラーキン法 要素マトリクス生成 宿題 プログラムの実行 データ構造 プログラムの構成 計算効率について McroAVS による可視化
24 Wndows 端末の使用法 小演習室 電源投入時 23 の 3 を選択する Mcro AVS と表示されています PuTTYによりun.ecc.u-too.ac.jp(ECCS28) に接続 実行 ファイルコピー McroAVS 立ち上げ Wndows ECCS28とファイルシステムを共有 文書ファイルは直接読むことができる
PuTTY 起動 (/4) / 25
PuTTY 起動 (2/4)SSH / 26
PuTTY 起動 (3/4) 27
PuTTY 起動 (4/4) / 28
29 McroAVS Wndows 専用 下記から評価版 (ヶ月有効 継続不可 ) をインストール可能 https://gt.cbernet.co.jp/v/dl/mavs_tral/mavs_tral.html
McroAVS co 起動 (/4) / 3
McroAVS co 起動 (2/4) / 3
32 McroAVS co 起動 (3/4) ドキュメント がECCS28における Documents に相当する
33 McroAVS co 起動 (4/4) / np ファイルが表示されていることを確認
34 UCDフォーマットについて (/4) Unstructured Cell Data 要素の種類キーワード 点 線 三角形 四角形 四面体 角錐 三角柱 六面体 二次要素 線 2 三角形 2 四角形 2 四面体 2 pt lne tr quad tet pr prsm he lne2 tr2 quad2 tet2 角錐 2 pr2 三角柱 2 六面体 2 prsm2 he2
35 UCDフォーマットについて (2/4) 書式の概要 ファイルの拡張子 データファイルの拡張子は.np np 書式 アスキーファイル 時系列 ( 複数ステップ ) に対応したものが標準フォーマット McroAVS Ver.6.( 現在は Ver..) ) まで使用していた旧フォーマット ( 単一ステップデータのための書式 ) も読み込むことは可能 実は一部これを使用している
36 UCDフォーマットについて (4/3) 書式の概要 ( コメント行 ) ( 要素のデータ成分数 )( 成分 の構成数 )( 成分 2 の構成数 ) ( 各成分の構成数 ) ( ステップ数 ) ( 要素データ成分 のラベル )( 単位 ) ( データの繰り返しタイプ ) ( 要素データ成分 2のラベル )( 単位 ) ( ステップ番号 ) ( コメント ) ( 全節点数 ) ( 全要素数 ) ( 節点番号 ) ( 座標 ) (Y 座標 ) (Z 座標 ) ( 節点番号 2) ( 座標 ) (Y 座標 ) (Z 座標 ) ( 要素番号 ) ( 材料番号 ) ( 要素の種類 ) ( 要素を構成する節点のつながり ) ( 要素番号 2) )( 材料番号 )( 要素の種類 )( 要素を構成する節点のつながり ) ( 各節点のデータ数 ) ( 各要素のデータ数 ) ( 節点のデータ成分数 )( 成分 の構成数 )( 成分 2 の構成数 ) ( 各成分の構成数 ) ( 節点データ成分 のラベル )( 単位 ) ( 節点データ成分 2のラベル )( 単位 ) ( 各要素データ成分のラベル )( 単位 ) ( 要素番号 ) ( 要素データ ) ( 要素データ 2) ( 要素番号 2) ( 要素データ ) ( 要素データ 2) ( ステップ番号 2) ( コメント ) ( 全節点数 ) ( 全要素数 ) ( 各節点データ成分のラベル )( 単位 ) ( 節点番号 ) ( 節点データ ) ( 節点データ 2) ( 節点番号 2) ( 節点データ ) ( 節点データ 2)
37 UCDフォーマットについて (4/4) 旧フォーマット ( 全節点数 ) ( 全要素数 ) ( 各節点のデータ数 ) ( 各要素のデータ数 ) ( モデルのデータ数 ) ( 節点番号 ) ( 座標 ) (Y 座標 ) (Z 座標 ) ( 節点番号 2) ( 座標 ) (Y 座標 ) (Z 座標 ) ( 要素番号 ) ( 材料番号 ) ( 要素の種類 ) ( 要素を構成する節点のつながり ) ( 要素番号 2) ( 材料番号 ) ( 要素の種類 ) ( 要素を構成する節点のつながり ) ( 要素のデータ成分数 ) ( 成分 の構成数 ) ( 成分 2の構成数 ) ( 各成分の構成数 ) ( 要素データ成分 のラベル )( 単位 ) ( 要素データ成分 2のラベル )( 単位 ) ( 各要素データ成分のラベル )( 単位 ) ( 要素番号 ) ( 要素データ ) ( 要素データ 2) ( 要素番号 2) ( 要素データ ) ( 要素データ 2) ( 節点のデータ成分数 ) ( 成分 の構成数 ) ( 成分 2 の構成数 ) ( 各成分の構成数 ) ( 節点データ成分 のラベル )( 単位 ) ( 節点データ成分 2 のラベル )( 単位 ) ( 各節点データ成分のラベル )( 単位 ) ( 節点番号 ) ( 節点データ) ( 節点データ2) ( 節点番号 2) ( 節点データ ) ( 節点データ 2)
38 UCDデータの例 # data 新フォーマット step 36 2... 2... 3 2... 4... 5... 6 2... 3.. 3. 32.. 3. 33 2.. 3. 34. 2. 3. 35. 2. 3. 36 2. 2. 3. he 2 5 4 4 3 2 he 2 3 6 5 2 5 4 3 he 4 5 8 7 3 4 7 6 he 2 2 24 23 29 3 33 32 he 22 23 26 25 3 32 35 34 2 he 23 24 27 26 32 33 36 35 Mesh. 2. 3. 4... 2. 36 2... 2... 旧フォーマット 3 2... 4... 5... 6 2... 3.. 3. 32.. 3. 33 2.. 3. 34. 2. 3. 35. 2. 3. 36 2. 2. 3. he 2 5 4 4 3 2 he 2 3 6 5 2 5 4 3 he 4 5 8 7 3 4 7 6 he 2 2 24 23 29 3 33 32 he 22 23 26 25 3 32 35 34 2 he 23 24 27 26 32 33 36 35 Mesh. 2. 3. 4... 2.
あああ ファイルを読み込む ファイル (F) データの読み込み (R) FEM3D-Part 39
このように 表示窓 一杯に広がってしまったら 物体の表示窓 から 物体 (O) 位置の指定 (T) を選択する FEM3D-Part 4
このような 位置の指定 ウィンドウが 表示される 最適な大きさ ズーム ( 縮小 ) な 最適な大きさ ズーム ( 縮小 ) な どのボタンを使用 どのボタンを使用 FEM3D-Part 4
様々な機能 : クリックして選択 表面コンター 変形図 42
メッシュ表示 クリック選択 物体の表示窓 から 物体(O) 表示方 FEM3D-Part 法の指定 (M) を選択する 出てきたウインドウの 線と滑らか表示 をクリックする 43
リジェンド表示 メインウィンドウのこのボタンを押すとリジェンドが表示され設定用のウィンドウも現れる 44
ファイル保存 ファイル (F) データの書き込み (W) で様々なフォーマットでの保存可能 動画も可能 45