有限要素法による 三次元定常熱伝導解析プログラム Fortran 編 202 年夏季集中講義中島研吾 並列計算プログラミング (66-2057) 先端計算機演習 (66-4009)
FEM3D 2 対象とする問題 : 三次元定常熱伝導 Z T T=0@Z=z ma Z T z 定常熱伝導 + 発熱 一様な熱伝導率 直方体 T z 一辺長さの立方体 ( 六面体 ) 要素 各方向にX Y Z 個 境界条件 Q z 0 X Y X Y T=0@Z=z ma 体積当たり発熱量は位置 ( メッシュの中心の座標 c c ) に依存 Q z QVOL C C
FEM3D 3 対象とする問題 : 三次元定常熱伝導 T T z T z Q z 0 原点から遠い部分が高温 体積当たり発熱量は位置 ( メッシュの中心の座標 ) に依存 Q z C C move
FEM3D 4 有限要素法の処理 支配方程式 ガラーキン法 : 弱形式 要素単位の積分 要素マトリクス生成 全体マトリクス生成 境界条件適用 連立一次方程式
FEM3D 5 有限要素法の処理 : プログラム 初期化 制御変数読み込み 座標読み込み 要素生成 (: 節点数 ICELTOT: 要素数 ) 配列初期化 ( 全体マトリクス 要素マトリクス ) 要素 全体マトリクスマッピング (IndeItem) マトリクス生成 要素単位の処理 (do cel= ICELTOT) 要素マトリクス計算 全体マトリクスへの重ね合わせ 境界条件の処理 連立一次方程式 共役勾配法 (CG)
FEM3D 6 三次元要素の定式化 三次元熱伝導方程式 ガラーキン法 要素マトリクス生成 プログラムの実行 データ構造 プログラムの構成
FEM3D 7 二次元への拡張 : 三角形要素 任意の形状を扱うことができる 特に一次要素は精度が悪く 一部の問題を除いてあまり使用されない
FEM3D 二次元への拡張 : 四角形要素 一次元要素と同じ形状関数を 軸に適用することによって 四角形要素の定式化は可能である 三角形と比較して特に低次要素の精度はよい しかしながら 各辺が座標軸に平行な長方形でなければならない 差分法と変わらない 4 3 このような形状を扱うことができない 2
FEM3D 9 アイソパラメトリック要素 (/3) 各要素を 自然座標系 () の正方形要素 [± ±] に変換する 4 3 4 + 3 2 - + - 2 各要素の全体座標系 (global coordnate)() における座標成分を 自然座標系における形状関数 [] ( 従属変数の内挿に使うのと同じ []) を使用して変換する場合 このような要素をアイソパラメトリック要素 (soparametrc element) という
アイソパラメトリック要素 (2/3) 2 3 4 各節点の座標 :( ) ( 2 2 ) ( 3 3 ) ( 4 4 ) 各節点における X 方向変位 :T T 2 T 3 T 4 T T 4 4 4 ) ( ) ( ) ( 4 3 - + + - 2 FEM3D 0
FEM3D アイソパラメトリック要素 (3/3) + - + - 高次の補間関数を使えば 曲線 曲面も扱うことが可能となる そういう意味で 自然座標系 と呼んでいる Sub-Parametrc Super-Parametrc
2D 自然座標系の形状関数 (/2) 自然座標系における正方形上の内挿多項式は下式で与えられる : 4 4 4 4 4 3 2 4 4 3 2 3 4 3 2 2 4 3 2 T T T T T T T T T T T T T T T T 4 3 2 T 各節点での条件より : 4 - + + - 3 2 FEM3D 2
2D 自然座標系の形状関数 (2/2) 元の式に代入して T について整理すると以下のようになる : 4 ) ( 4 ) ( 4 ) ( 4 ) ( 4 3 2 形状関数 は以下のようになる : 双一次 (b-lnear) 要素とも呼ばれる 各節点における の値を計算してみよ - + + - 4 3 2 FEM3D 3 4 4 3 3 2 2 T T T T T
FEM3D 4 三次元への拡張 四面体要素 : 二次元における三角形要素 任意の形状を扱うことができる 特に一次要素は精度が悪く 一部の問題を除いてあまり使用されない 高次の四面体要素は広く使用されている 本講義では低次六面体要素 ( アイソパラメトリック要素 ) を使用する tr-lnear 自由度 : 温度 @ 各節点上
3D 自然座標系の形状関数 ) ( ) ( ) ( ) ( 4 3 2 2 3 4 5 6 7 ) ( ) ( ) ( ) ( 7 6 5 w w v v u u ) ( ) ( ) ( FEM3D 5
FEM3D 6 三次元要素の定式化 三次元熱伝導方程式 ガラーキン法 要素マトリクス生成 プログラムの実行 データ構造 プログラムの構成
ガラーキン法の適用 (/3) 以下のような三次元熱伝導方程式を考慮する ( 熱伝導率一定 ): } { T 要素内の温度分布 ( マトリクス形式 ) 節点における温度を としてある ガラーキン法に従い 重み関数を [] とすると 各要素において以下の積分方程式が得られる : 0 2 2 2 2 2 2 dv Q z T T T T V 0 2 2 2 2 2 2 Q z T T T FEM3D 7
ガラーキン法の適用 (2/3) 三次元のグリーンの定理 前式の 2 階微分の部分に適用 ( 表面積分省略 ): これに以下を代入する : } { } { } { } { z z T T T T V V S dv z B z A B A B A ds n B A dv z B B B A 2 2 2 2 2 2 dv T T T dv T T T V z T z T T zz T V FEM3D
FEM3D 9 ガラーキン法の適用 (3/3) 体積あたり発熱量の項を加えて次式が得られる : Q T T T z z dv Q dv 0 V この式を弱形式 (weak form) と呼ぶ 元の微分方程式では 2 階の微分が含まれていたが 上式では グリーンの定理によって 階微分に低減されている 弱形式によって近似関数 ( 形状関数 内挿関数 ) に対する要求が弱くなっている : すなわち線形関数で2 階微分の効果を記述できる 項が増えただけで 一次元と同じ V
境界条件を考慮した弱形式 : 各要素 dv dv dv k V z T z V T V T e ) ( ) ( ) ( ) ( e e e f k dv Q f V T e ) ( FEM3D 20
FEM3D 2 要素マトリクス : 行列 j j k j 7 5 6 4 3 2
FEM3D 22 要素マトリクス :k j j k j 7 5 6 4 3 2 dv T ( e) k T V V T z z dv V dv k j V j j z j z dv
あとは積分を求めれば良い 自然座標系 () 上で定義 ガウス積分公式が使える ( 三次元 ) しかし 微分が k j W W W ) ( k j LM: 方向の積分点数 : 積分点での重み係数 : 積分点の座標値 k k j k j M j L f W W W d d d f I ) ( ) ( FEM3D 23
自然座標系における偏微分 (/4) 偏微分の公式より以下のようになる : z z z z z z ) ( ) ( ) ( z は定義より簡単に求められるがを実際の計算で使用する FEM3D 24
自然座標系における偏微分 (2/4) マトリックス表示すると : z J z z z z J : ヤコビのマトリクス (Jacob matr Jacoban) FEM3D 25
自然座標系における偏微分 (3/4) の定義より簡単に求められる z z z J J J z z z J J J z z z J J J 33 32 3 23 22 2 3 2 FEM3D 26
自然座標系における偏微分 (4/4) 従って下記のように偏微分を計算できる ヤコビアン (3 3 行列 ) の逆行列を求める J z z z z FEM3D 27
要素単位での積分 V j j j V z j z j j j dv z z dv k 2 3 4 5 6 7 j k j FEM3D 2
自然座標系での積分 d d d J z z dddz z z dv z z j j j j j j V j j j det 2 3 4 5 6 7 j k j FEM3D 29
ガウスの積分公式 2 3 4 5 6 7 k k j k j M j L f W W W d d d f I ) ( ) ( d d d J z z j j j det FEM3D 30 j k j
FEM3D 3 残りの手順 ここまでで 要素ごとの積分が可能となる あとは : 全体マトリクスへの重ね合わせ 境界条件処理 連立一次方程式を解く 詳細はプログラムの内容を解説しながら説明する
要素 全体マトリクス重ね合わせ 2 4 5 6 3 2 2 3 4 2 2 4 3 () 4 () 3 () 2 () () 4 () 3 () 2 () () 44 () 43 () 42 () 4 () 34 () 33 () 32 () 3 () 24 () 23 () 22 () 2 () 4 () 3 () 2 () () () () } { } ]{ [ f f f f k k k k k k k k k k k k k k k k f k (2) 4 (2) 3 (2) 2 (2) (2) 4 (2) 3 (2) 2 (2) (2) 44 (2) 43 (2) 42 (2) 4 (2) 34 (2) 33 (2) 32 (2) 3 (2) 24 (2) 23 (2) 22 (2) 2 (2) 4 (2) 3 (2) 2 (2) (2) (2) (2) } { } ]{ [ f f f f k k k k k k k k k k k k k k k k f k 6 5 4 3 2 6 5 4 3 2 6 5 4 3 2 } { } ]{ [ B B B B B B D X X X X D X X X X D X X X X X X D X X X X D X X X X D F K FEM3D 32
FEM3D 33 三次元要素の定式化 三次元弾性力学方程式 ガラーキン法 要素マトリクス生成 演習 プログラムの実行 データ構造 プログラムの構成
FEM3D 34 演習 ガウスの積分公式を使用して以下の四角形の面積を求めよ ( プログラムを作って 計算機で計算してください ) 4 V 3 2 :(.0.0) 2:(4.0 2.0) 3:(3.0 5.0) 4:(2.0 4.0) I V dv det J dd
FEM3D 35 ヒント (/2) 座標値によってヤコビアン ( ヤコビの行列 ) を計算 ガウスの積分公式 (n=2) に代入する I m W W j f ( j ) f ( ) dd n j mplct REAL* (A-HO-Z) real* W(2) real* POI(2) W()=.0d0 W(2)=.0d0 POI()= -0.5773502692d0 POI(2)= +0.5773502692d0 SUM= 0.d0 do jp= 2 do p= 2 FC = F(POI(p)POI(jp)) SUM= SUM + W (p)*w (jp)*fc
ヒント (2/2) J J det 4 ) ( 4 ) ( 4 ) ( 4 ) ( 4 3 2 4 4 4 4 4 4 4 4 FEM3D 36
FEM3D 37 三次元要素の定式化 三次元弾性力学方程式 ガラーキン法 要素マトリクス生成 プログラムの実行 データ構造 プログラムの構成
FEM3D 3 対象とする問題 : 三次元定常熱伝導 Z T T=0@Z=z ma Z T z 定常熱伝導 + 発熱 一様な熱伝導率 直方体 T z 一辺長さの立方体 ( 六面体 ) 要素 各方向にX Y Z 個 境界条件 Q z 0 X Y X Y T=0@Z=z ma 体積当たり発熱量は位置 ( メッシュの中心の座標 c c ) に依存 Q z QVOL C C
FEM3D 39 ファイルコピー インストール 三次元熱伝導解析コード >$ cd <$E-TOP> >$ cp /home03/skengon/documents/class_eps/f/fem3d.tar. >$ cp /home03/skengon/documents/class_eps/c/fem3d.tar. >$ tar vf fem3d.tar >$ cd fem3d >$ ls run src インストール >$ cd <$E-TOP>/fem3d/src >$ make >$ ls../run/sol sol メッシュジェネレータインストール >$ cd <$E-TOP>/fem3d/run >$ g95 O3 mgcube.f o mgcube
FEM3D 40 計算の流れメッシュ生成 計算 ファイル名称固定 mgcube メッシュジェネレータ cube.0 メッシュファイル sol FEM 本体 IPUT.DAT 制御データ test.np 可視化用出力
FEM3D 4 メッシュ生成 Z >$ cd <$E-TOP>/fem3d/run >$./mgcube T=0@Z=z ma X Y Z 各辺長さを訊いてくる 202020 このように入れてみる Z >$ ls cube.0 生成を確認 cube.0 X Y X Y
FEM3D 42 制御ファイル :IPUT.DAT IPUT.DAT cube.0 fname 2000 ITER.0.0 COD QVOL.0e-0 RESID fname: ITER: COD: QVOL: RESID: メッシュファイル名反復回数上限熱伝導率体積当たり発熱量係数反復法の収束判定値 T Q T z z QVOL C C T z Q z 0
FEM3D 43 >$ cd <$E-TOP>/fem3d/run >$./sol 実行 >$ ls test.np 生成を確認 test.np
FEM3D 44 ParaVew http://www.paravew.org/ ファイルを開く 図の表示 イメージファイルの保存 http://nkl.cc.utoko.ac.jp/class/howtouseparavew.pdf
FEM3D 45 UCD Format (/3) Unstructured Cell Data 要素の種類点線三角形四角形四面体角錐三角柱六面体二次要素線 2 三角形 2 四角形 2 四面体 2 角錐 2 三角柱 2 六面体 2 キーワード pt lne tr quad tet pr prsm he lne2 tr2 quad2 tet2 pr2 prsm2 he2
FEM3D 46 UCD Format (2/3) Orgnall for AVS mcroavs Etenson of the UCD fle s np There are two tpes of formats. Onl old tpe can be read b ParaVew.
FEM3D 47 UCD Format (3/3): Old Format ( 全節点数 ) ( 全要素数 ) ( 各節点のデータ数 ) ( 各要素のデータ数 ) ( モデルのデータ数 ) ( 節点番号 ) (X 座標 ) (Y 座標 ) (Z 座標 ) ( 節点番号 2) (X 座標 ) (Y 座標 ) (Z 座標 ) ( 要素番号 ) ( 材料番号 ) ( 要素の種類 ) ( 要素を構成する節点のつながり ) ( 要素番号 2) ( 材料番号 ) ( 要素の種類 ) ( 要素を構成する節点のつながり ) ( 要素のデータ成分数 ) ( 成分 の構成数 ) ( 成分 2 の構成数 ) ( 各成分の構成数 ) ( 要素データ成分 のラベル )( 単位 ) ( 要素データ成分 2 のラベル )( 単位 ) ( 各要素データ成分のラベル )( 単位 ) ( 要素番号 ) ( 要素データ ) ( 要素データ 2) ( 要素番号 2) ( 要素データ ) ( 要素データ 2) ( 節点のデータ成分数 ) ( 成分 の構成数 ) ( 成分 2 の構成数 ) ( 各成分の構成数 ) ( 節点データ成分 のラベル )( 単位 ) ( 節点データ成分 2 のラベル )( 単位 ) ( 各節点データ成分のラベル )( 単位 ) ( 節点番号 ) ( 節点データ ) ( 節点データ 2) ( 節点番号 2) ( 節点データ ) ( 節点データ 2)
FEM3D 4 三次元要素の定式化 三次元弾性力学方程式 ガラーキン法 要素マトリクス生成 プログラムの実行 データ構造 プログラムの構成
FEM3D 49 メッシュファイル構成 :cube.0 番号は から始まっている 節点データ 節点数 節点番号 座標 要素データ 要素数 要素タイプ 要素番号 材料番号 コネクティビティ 節点グループデータ グループ数 グループ内節点数 グループ名 グループ内節点
FEM3D 50 メッシュ生成コード :mgcube.f(/5) mplct REAL* (A-HO-Z) real(knd=) dmenson(::) allocatable :: X Y real(knd=) dmenson(::) allocatable :: X0 Y0 character(len=0) :: GRIDFILE HHH nteger dmenson(::) allocatable :: IW!C!C +-------+!C IIT.!C +-------+!C=== wrte (**) 'XYZ' read (**) XYZ XP= X + YP= Y + ZP= Z + DX=.d0 IODTOT= XP*YP*ZP ICELTOT= X *Y *Z IBODTOT= XP*YP allocate (IW(IODTOT4)) IW= 0 X 方向節点数 Y 方向節点数 Z 方向節点数 総節点数総要素数 XY 平面の節点数
FEM3D 5 メッシュ生成コード :mgcube.f(2/5) cou= 0!C=== b = do k= ZP do j= YP = cou= cou + = (k-)*ibodtot + (j-)*xp + IW(coub)= cou= 0 b = 2 do k= ZP j= do = XP cou= cou + = (k-)*ibodtot + (j-)*xp + IW(coub)= cou= 0 b = 3 k= do j= YP do = XP cou= cou + = (k-)*ibodtot + (j-)*xp + IW(coub)= cou= 0 b = 4 k= ZP do j= YP do = XP cou= cou + = (k-)*ibodtot + (j-)*xp + IW(coub)= X=Xmn の節点を IW(b) に格納 (b=yp*zp) X Z T=0@Z=z ma Y X Z Y
FEM3D 52 メッシュ生成コード :mgcube.f(2/5) cou= 0 b = do k= ZP do j= YP = cou= cou + = (k-)*ibodtot + (j-)*xp + IW(coub)=!C=== cou= 0 b = 2 do k= ZP j= do = XP cou= cou + = (k-)*ibodtot + (j-)*xp + IW(coub)= cou= 0 b = 3 k= do j= YP do = XP cou= cou + = (k-)*ibodtot + (j-)*xp + IW(coub)= cou= 0 b = 4 k= ZP do j= YP do = XP cou= cou + = (k-)*ibodtot + (j-)*xp + IW(coub)= Y=Ymn の節点を IW(b2) に格納 (b=xp*zp) X Z T=0@Z=z ma Y X Z Y
FEM3D 53 メッシュ生成コード :mgcube.f(2/5)!c=== cou= 0 b = do k= ZP do j= YP = cou= cou + = (k-)*ibodtot + (j-)*xp + IW(coub)= cou= 0 b = 2 do k= ZP j= do = XP cou= cou + = (k-)*ibodtot + (j-)*xp + IW(coub)= cou= 0 b = 3 k= do j= YP do = XP cou= cou + = (k-)*ibodtot + (j-)*xp + IW(coub)= cou= 0 b = 4 k= ZP do j= YP do = XP cou= cou + = (k-)*ibodtot + (j-)*xp + IW(coub)= X Z T=0@Z=z ma Y X Z Z=Zmn の節点を IW(b3) に格納 (b=xp*yp) Z=Zma の節点を IW(b4) に格納 (b=xp*yp) Y
FEM3D 54 メッシュ生成コード :mgcube.f(3/5)!c!c +-------------+!C GeoFEM data!c +-------------+!C=== = 0 wrte (**) 'GeoFEM grdfle name?' GRIDFILE= 'cube.0' open (2 fle= GRIDFILE status='unknown'form='formatted') wrte(2'(00)') IODTOT cou= 0 do k= ZP do j= YP do = XP XX= dfloat(-)*dx YY= dfloat(j-)*dx ZZ= dfloat(k-)*dx cou= cou + wrte (2'(03(pe6.6))') cou XX YY ZZ wrte(2'(0)') ICELTOT IELMTYPL= 36 wrte(2'(00)') (IELMTYPL =ICELTOT) 節点数節点番号 座標 要素数要素タイプ ( 使わないデータ )
FEM3D 55 cube.0: 節点データ 要素数 要素タイプ (X=Y=Z=4) 25 =5*5*5 0.000000E+00 0.000000E+00 0.000000E+00 2.000000E+00 0.000000E+00 0.000000E+00 3 2.000000E+00 0.000000E+00 0.000000E+00 4 3.000000E+00 0.000000E+00 0.000000E+00 5 4.000000E+00 0.000000E+00 0.000000E+00 6 0.000000E+00.000000E+00 0.000000E+00 7.000000E+00.000000E+00 0.000000E+00 2.000000E+00.000000E+00 0.000000E+00 9 3.000000E+00.000000E+00 0.000000E+00 ( 途中省略 ) 2 0.000000E+00 4.000000E+00 4.000000E+00 22.000000E+00 4.000000E+00 4.000000E+00 23 2.000000E+00 4.000000E+00 4.000000E+00 24 3.000000E+00 4.000000E+00 4.000000E+00 25 4.000000E+00 4.000000E+00 4.000000E+00 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 56 メッシュ生成コード :mgcube.f(4/5) cou= 0 mat= do k= Z do j= Y do = X cou= cou + n = (k-)*ibodtot + (j-)*xp + n2 = n + n3 = n2 + XP n4 = n3 - n5 = n + IBODTOT n6 = n2 + IBODTOT n7 = n3 + IBODTOT n = n4 + IBODTOT wrte (2'(00)') cou mat n n2 n3 n4 & & n5 n6 n7 n mat: 材料番号 (=) 7 5 6 4 3 2
FEM3D 57 cube.0: 要素データ 2 7 6 26 27 32 3 2 2 3 7 27 2 33 32 3 3 4 9 2 29 34 33 4 4 5 0 9 29 30 35 34 5 6 7 2 3 32 37 36 6 7 3 2 32 33 3 37 7 9 4 3 33 34 39 3 9 0 5 4 34 35 40 39 9 2 7 6 36 37 42 4 0 2 3 7 37 3 43 42 3 4 9 3 39 44 43 2 4 5 20 9 39 40 45 44 3 6 7 22 2 4 42 47 46 ( 途中省略 ) 42 62 63 6 67 7 93 92 43 63 64 69 6 9 94 93 44 64 65 70 69 9 90 95 94 45 66 67 72 7 9 92 97 96 46 67 6 73 72 92 93 9 97 47 6 69 74 73 93 94 99 9 4 69 70 75 74 94 95 00 99 49 76 77 2 0 02 07 06 50 77 7 3 2 02 03 0 07 5 7 79 4 3 03 04 09 0 52 79 0 5 4 04 05 0 09 53 2 7 6 06 07 2 54 2 3 7 07 0 3 2 55 3 4 9 0 09 4 3 56 4 5 90 9 09 0 5 4 57 6 7 92 9 2 7 6 5 7 93 92 2 3 7 59 9 94 93 3 4 9 60 9 90 95 94 4 5 20 9 6 9 92 97 96 6 7 22 2 62 92 93 9 97 7 23 22 63 93 94 99 9 9 24 23 64 94 95 00 99 9 20 25 24
FEM3D 5 メッシュ生成コード :mgcube.f(5/5) IGTOT= 4 IBT= YP*ZP IBT2= XP*ZP + IBT IBT3= XP*YP + IBT2 IBT4= XP*YP + IBT3 wrte (2'(00)') IGTOT wrte (2'(00)') IBT IBT2 IBT3 IBT4 HHH= 'Xmn' wrte (2'(a0)') HHH wrte (2'(00)') (IW() =YP*ZP) HHH= 'Ymn' wrte (2'(a0)') HHH wrte (2'(00)') (IW(2) =XP*ZP) HHH= 'Zmn' wrte (2'(a0)') HHH wrte (2'(00)') (IW(3) =XP*YP) HHH= 'Zma' wrte (2'(a0)') HHH wrte (2'(00)') (IW(4) =XP*YP) IGTOT IBT グループ総数 (XmnYmnZmnZma) 累積数 ( 以下略 )!C=== stop end deallocate (IW) close (2)
FEM3D 59 cube.0: 節点グループデータ Xmn Ymn Zmn Zma 4 25 50 75 00 6 6 2 26 3 36 4 46 5 56 6 66 7 76 6 9 96 0 06 6 2 2 3 4 5 26 27 2 29 30 5 52 53 54 55 76 77 7 79 0 0 02 03 04 05 2 3 4 5 6 7 9 0 2 3 4 5 6 7 9 20 2 22 23 24 25 0 02 03 04 05 06 07 0 09 0 2 3 4 5 6 7 9 20 2 22 23 24 25 ( 以下使用せず )
FEM3D 60 メッシュ生成 実は技術的には大きな課題 複雑形状 大規模メッシュ 並列化が難しい 市販のメッシュ生成アプリケーション FEMAP CAD データとのインタフェース move
FEM3D 6 三次元要素の定式化 三次元弾性力学方程式 ガラーキン法 要素マトリクス生成 プログラムの実行 データ構造 プログラムの構成
FEM3D 62 有限要素法の処理 : プログラム 初期化 制御変数読み込み 座標読み込み 要素生成 (: 節点数 E: 要素数 ) 配列初期化 ( 全体マトリクス 要素マトリクス ) 要素 全体マトリクスマッピング (IndeItem) マトリクス生成 要素単位の処理 (do cel= E) 要素マトリクス計算 全体マトリクスへの重ね合わせ 境界条件の処理 連立一次方程式 共役勾配法 (CG)
FEM3D 63 test メインプログラム nput_cntl 制御データ入力 nput_grd メッシュファイル入力 fnd_node 節点探索 mat_con0 行列コネクティビティ生成 msort ソート mat_con 行列コネクティビティ生成 mat_ass_man 係数行列生成 jacob ヤコビアン計算 mat_ass_bc 境界条件処理 三次元熱伝導解 析コード heat3d の構成 solve 線形ソルバー制御 output_ucd 可視化処理 cg CG 法計算
FEM3D 64 全体処理 program heat3d use solver use pfem_utl mplct REAL*(A-HO-Z) call IPUT_CTL call IPUT_GRID call MAT_CO0 call MAT_CO call MAT_ASS_MAI call MAT_ASS_BC call SOLVE call OUTPUT_UCD end program heat3d
FEM3D 65 Global 変数表 :pfem_utl.f(/3) 変数名種別サイズ I/O 内容 fname C (0) I メッシュファイル名 P I I 節点数 ICELTOT I I 要素数 ODGRPtot I I 節点グループ数 XYZ R (3) I 節点座標 ICELOD I (ICELTOT) I 要素コネクティビティ ODGRP_IDEX I (0:ODGRPtot) I 各節点グループに含まれる節点数 ( 累積 ) ODGRP_ITEM I (ODGRP_IDEX(ODG RPTOT) I 節点グループに含まれる節点 ODGRP_AME C0 (ODGRPTOT) I 節点グループ名 LU I O 各節点非対角成分数 PLU I O 非対角成分総数 D R () O 全体行列 : 対角ブロック BX R () O 右辺ベクトル 未知数ベクトル
FEM3D 66 Global 変数表 :pfem_utl.f(2/3) 変数名種別サイズ I/O 内容 AMAT R (PLU) O 全体行列 : 非零非対角成分 nde I (0:) O 全体行列 : 非零非対角成分数 tem I (PLU) O 全体行列 : 非零非対角成分 ( 列番号 ) ILU I () O 各節点の非零非対角成分数 IALU I (LU) O 各節点の非零非対角成分数 ( 列番号 ) IWKX I (2) O ワーク用配列 ITER ITERactual I I 反復回数の上限 実際の反復回数 RESID R I 打ち切り誤差 (.e-に設定) pfemiarra I (00) O 諸定数 ( 整数 ) pfemrarra R (00) O 諸定数 ( 実数 )
FEM3D 67 Global 変数表 :pfem_utl.f(3/3) 変数名種別サイズ I/O 内容 Oth R I =0.25 PQ PE PT R (22) O 各ガウス積分点における POS WEI R (22) O 各ガウス積分点の座標 重み係数 COL COL2 I (00) O ソート用ワーク配列 SHAPE R (222) O 各ガウス積分点における形状関数 (=~) PX PY PZ R (222) O 各ガウス積分点における ~ ~ z DETJ R (222) O 各ガウス積分点におけるヤコビアン行列式 COD QVOL R I 熱伝導率 体積当たり発熱量係数
FEM3D 6 制御ファイル入力 :IPUT_CTL subroutne IPUT_CTL use pfem_utl mplct REAL* (A-HO-Z) open (fle= 'IPUT.DAT' status='unknown') read ('(a0)') fname read (*) ITER read (*) COD QVOL read (*) RESID close () pfemiarra()= ITER pfemrarra()= RESID return end IPUT.DAT cube.0 fname 2000 ITER.0.0 COD QVOL.0e-0 RESID
FEM3D 69 メッシュ入力 :IPUT_GRID(/2)!C***!C*** IPUT_GRID!C***!C subroutne IPUT_GRID use pfem_utl mplct REAL* (A-HO-Z) open ( fle= fname status= 'unknown' form= 'formatted')!c!c-- ODE read (*) P= allocate (XYZ(3)) XYZ= 0.d0 do = read (*) (XYZ(kk)kk=3)
FEM3D 70 allocate deallocate 関数 (C 言語 ) #nclude <stdo.h> #nclude <stdlb.h> vod* allocate_vector(nt szent m) { vod *a; f ( ( a=(vod * )malloc( m * sze ) ) == ULL ) { fprntf(stdout"error:memor does not enough! n vector n"); et(); } return a; } vod deallocate_vector(vod *a) { free( a ); } vod** allocate_matr(nt szent mnt n) { vod **aa; nt ; f ( ( aa=(vod ** )malloc( m * szeof(vod*) ) ) == ULL ) { fprntf(stdout"error:memor does not enough! aa n matr n"); et(); } f ( ( aa[0]=(vod * )malloc( m * n * sze ) ) == ULL ) { fprntf(stdout"error:memor does not enough! n matr n"); et(); } for(=;<m;++) aa[]=(char*)aa[-]+sze*n; return aa; } vod deallocate_matr(vod **aa) { free( aa ); } allocate を FORTRA 並みに簡単にやるための関数
FEM3D 7 メッシュ入力 :IPUT_GRID(2/2)!C!C-- ELEMET read (*) ICELTOT allocate (ICELOD(ICELTOT)) read ('(00)') (TYPE = ICELTOT) do cel= ICELTOT read ('(0025)') IMAT (ICELOD(celk) k=)!c!c-- ODE grp. nfo. read ('(00)') ODGRPtot allocate (ODGRP_IDEX(0:ODGRPtot)ODGRP_AME(ODGRPtot)) ODGRP_IDEX= 0 read ('(00)') (ODGRP_IDEX() = ODGRPtot) nn= ODGRP_IDEX(ODGRPtot) allocate (ODGRP_ITEM(nn)) do k= ODGRPtot S= ODGRP_IDEX(k-) + E= ODGRP_IDEX(k ) read ('(a0)') ODGRP_AME(k) nn= E - S + f (nn.ne.0) then read ('(00)') (ODGRP_ITEM(kk)kk=S E) endf close () return end
FEM3D 72 test メインプログラム nput_cntl 制御データ入力 nput_grd メッシュファイル入力 fnd_node 節点探索 mat_con0 行列コネクティビティ生成 msort ソート mat_con 行列コネクティビティ生成 mat_ass_man 係数行列生成 jacob ヤコビアン計算 mat_ass_bc 境界条件処理 三次元熱伝導解 析コード heat3d の構成 solve 線形ソルバー制御 output_ucd 可視化処理 cg CG 法計算
FEM3D 73 Global 変数表 :pfem_utl.h/f(/3) 変数名種別サイズ I/O 内容 fname C (0) I メッシュファイル名 P I I 節点数 ICELTOT I I 要素数 ODGRPtot I I 節点グループ数 XYZ R (3) I 節点座標 ICELOD I (ICELTOT) I 要素コネクティビティ ODGRP_IDEX I (0:ODGRPtot) I 各節点グループに含まれる節点数 ( 累積 ) ODGRP_ITEM I (ODGRP_IDEX(ODG RPTOT) I 節点グループに含まれる節点 ODGRP_AME C0 (ODGRPTOT) I 節点グループ名 LU I O 各節点非対角成分数 PLU I O 非対角成分総数 D R () O 全体行列 : 対角ブロック BX R () O 右辺ベクトル 未知数ベクトル
FEM3D 74 Global 変数表 :pfem_utl.h/f(2/3) 変数名種別サイズ I/O 内容 AMAT R (PLU) O 全体行列 : 非零非対角成分 nde I (0:) O 全体行列 : 非零非対角成分数 tem I (PLU) O 全体行列 : 非零非対角成分 ( 列番号 ) ILU I () O 各節点の非零非対角成分数 IALU I (LU) O 各節点の非零非対角成分数 ( 列番号 ) IWKX I (2) O ワーク用配列 ITER ITERactual I I 反復回数の上限 実際の反復回数 RESID R I 打ち切り誤差 (.e-に設定) pfemiarra I (00) O 諸定数 ( 整数 ) pfemrarra R (00) O 諸定数 ( 実数 )
FEM3D 75 Global 変数表 :pfem_utl.h/f(3/3) 変数名種別サイズ I/O 内容 Oth R I =0.25 PQ PE PT R (22) O 各ガウス積分点における POS WEI R (22) O 各ガウス積分点の座標 重み係数 COL COL2 I (00) O ソート用ワーク配列 SHAPE R (222) O 各ガウス積分点における形状関数 (=~) PX PY PZ R (222) O 各ガウス積分点における ~ ~ z DETJ R (222) O 各ガウス積分点におけるヤコビアン行列式 COD QVOL R I 熱伝導率 体積当たり発熱量係数
FEM3D 76 マトリクス生成まで 一次元のときは ndetem に関連した情報を簡単に作ることができた 非ゼロ非対角成分の数は2 番号が自分に対して :+と- 三次元の場合はもっと複雑 非ゼロ非対角成分の数は7~26( 現在の形状 ) 実際はもっと複雑 前以て 非ゼロ非対角成分数はわからない move
FEM3D 77 マトリクス生成まで 一次元のときは ndetem に関連した情報を簡単に作ることができた 非ゼロ非対角成分の数は2 番号が自分に対して :+と- 三次元の場合はもっと複雑 非ゼロ非対角成分の数は7~26( 現在の形状 ) 実際はもっと複雑 前以て 非ゼロ非対角成分の数はわからない ILU()IALU(LU) を使って非ゼロ非対角成分数を予備的に勘定する
FEM3D 7 全体処理 program heat3d use solver use pfem_utl mplct REAL*(A-HO-Z) call IPUT_CTL call IPUT_GRID call MAT_CO0 call MAT_CO call MAT_ASS_MAI call MAT_ASS_BC call SOLVE call OUTPUT_UCD 3 4 5 6 7 9 9 0 2 4 5 6 5 6 7 2 3 2 3 4 end program heat3d MAT_CO0: IU IALU 生成 MAT_CO: nde tem 生成 とりあえず から始まる節点番号を記憶
FEM3D 79 MAT_CO0: 全体構成 do cel= ICELTOT 節点相互の関係から ILUIALU を生成 (FID_ODE) 7 5 6 4 3 3 4 5 6 7 9 9 0 2 2 4 5 6 5 6 7 2 3 2 3 4
FEM3D 0 行列コネクティビティ生成 : MAT_CO0(/4)!C!C***!C*** MAT_CO0!C***!C subroutne MAT_CO0 use pfem_utl mplct REAL* (A-HO-Z) LU= 26 allocate (ILU() IALU(LU)) ILU= 0 IALU= 0 LU: 各節点における非ゼロ非対角成分の最大数 ( 接続する節点数 ) 今の問題の場合はわかっているので このようにできる 不明の場合の実装 : レポート課題
FEM3D 行列コネクティビティ生成 : MAT_CO0(/4)!C!C***!C*** MAT_CO0!C***!C subroutne MAT_CO0 use pfem_utl mplct REAL* (A-HO-Z) LU= 26 allocate (ILU() IALU(LU)) ILU= 0 IALU= 0 変数名サイズ内容 ILU IALU () (LU) 各節点の非零非対角成分数 各節点の非零非対角成分 ( 列番号 )
FEM3D 2 行列コネクティビティ生成 : MAT_CO0(2/4): から始まる番号 do cel= ICELTOT n= ICELOD(cel) n2= ICELOD(cel2) n3= ICELOD(cel3) n4= ICELOD(cel4) n5= ICELOD(cel5) n6= ICELOD(cel6) n7= ICELOD(cel7) n= ICELOD(cel) call FID_TS_ODE (nn2) call FID_TS_ODE (nn3) call FID_TS_ODE (nn4) call FID_TS_ODE (nn5) call FID_TS_ODE (nn6) call FID_TS_ODE (nn7) call FID_TS_ODE (nn) 7 5 6 4 3 2 call FID_TS_ODE (n2n) call FID_TS_ODE (n2n3) call FID_TS_ODE (n2n4) call FID_TS_ODE (n2n5) call FID_TS_ODE (n2n6) call FID_TS_ODE (n2n7) call FID_TS_ODE (n2n) call FID_TS_ODE (n3n) call FID_TS_ODE (n3n2) call FID_TS_ODE (n3n4) call FID_TS_ODE (n3n5) call FID_TS_ODE (n3n6) call FID_TS_ODE (n3n7) call FID_TS_ODE (n3n)
FEM3D 3 節点探索 :FID_TS_ODE ILUIAU 探索 : 一次元ではこの部分は手動!C!C***!C*** FID_TS_ODE!C***!C subroutne FID_TS_ODE (pp2) do kk= ILU(p) f (p2.eq.ialu(pkk)) return cou= ILU(p) + IALU(pcou)= p2 ILU(p )= cou return end subroutne FID_TS_ODE 変数名サイズ内容 ILU IALU () (LU) 各節点の非零非対角成分数 各節点の非零非対角成分 ( 列番号 )
FEM3D 4 節点探索 :FID_TS_ODE 一次元ではこの部分は手動!C!C***!C*** FID_TS_ODE!C***!C subroutne FID_TS_ODE (pp2) do kk= ILU(p) f (p2.eq.ialu(pkk)) return cou= ILU(p) + IALU(pcou)= p2 ILU(p )= cou return 既に IALU に含まれている場合は 次のペアへ end subroutne FID_TS_ODE 3 4 5 6 7 9 9 0 2 4 5 6 5 6 7 2 3 2 3 4
FEM3D 5 節点探索 :FID_TS_ODE 一次元ではこの部分は手動!C!C***!C*** FID_TS_ODE!C***!C subroutne FID_TS_ODE (pp2) do kk= ILU(p) f (p2.eq.ialu(pkk)) return cou= ILU(p) + IALU(pcou)= p2 ILU(p )= cou return end subroutne FID_TS_ODE IALU に含まれていない場合は ILU に を加えて IALU に格納 3 4 5 6 7 9 9 0 2 4 5 6 5 6 7 2 3 2 3 4
FEM3D 6 行列コネクティビティ生成 : MAT_CO0(3/4) call FID_TS_ODE (n4n) call FID_TS_ODE (n4n2) call FID_TS_ODE (n4n3) call FID_TS_ODE (n4n5) call FID_TS_ODE (n4n6) call FID_TS_ODE (n4n7) call FID_TS_ODE (n4n) call FID_TS_ODE (n5n) call FID_TS_ODE (n5n2) call FID_TS_ODE (n5n3) call FID_TS_ODE (n5n4) call FID_TS_ODE (n5n6) call FID_TS_ODE (n5n7) call FID_TS_ODE (n5n) call FID_TS_ODE (n6n) call FID_TS_ODE (n6n2) call FID_TS_ODE (n6n3) call FID_TS_ODE (n6n4) call FID_TS_ODE (n6n5) call FID_TS_ODE (n6n7) call FID_TS_ODE (n6n) 7 5 6 4 3 2 call FID_TS_ODE (n7n) call FID_TS_ODE (n7n2) call FID_TS_ODE (n7n3) call FID_TS_ODE (n7n4) call FID_TS_ODE (n7n5) call FID_TS_ODE (n7n6) call FID_TS_ODE (n7n)
FEM3D 7 行列コネクティビティ生成 : MAT_CO0(4/4) call FID_TS_ODE (nn) call FID_TS_ODE (nn2) call FID_TS_ODE (nn3) call FID_TS_ODE (nn4) call FID_TS_ODE (nn5) call FID_TS_ODE (nn6) call FID_TS_ODE (nn7) do n= = ILU(n) do k= COL(k)= IALU(nk) call msort (COL COL2 ) do k= - IALU(n-k+)= COL(COL2(k)) 各節点において IALU(k) が小さい番号から大きい番号に並ぶようにソート ( 単純なバブルソート ) せいぜい 00 程度のものをソートする
FEM3D CRS 形式への変換 :MAT_CO!C!C***!C*** MAT_CO!C***!C subroutne MAT_CO use pfem_utl mplct REAL* (A-HO-Z) allocate (nde(0:)) nde= 0 do = nde()= nde(-) + ILU() PLU= nde() allocate (tem(plu)) do = do k= ILU() kk = k + nde(-) tem(kk)= IALU(k) deallocate (ILU IALU) end subroutne MAT_CO C nde[ ] nde[0] 0 FORTRA nde(0) 0 k 0 ILU[ k] nde( ) ILU( k) k
FEM3D 9 CRS 形式への変換 :MAT_CO!C!C***!C*** MAT_CO!C***!C subroutne MAT_CO use pfem_utl mplct REAL* (A-HO-Z) allocate (nde(0:)) nde= 0 do = nde()= nde(-) + ILU() PLU= nde() allocate (tem(plu)) do = do k= ILU() kk = k + nde(-) tem(kk)= IALU(k) PLU=nde() tem のサイズ非ゼロ非対角成分総数 deallocate (ILU IALU) end subroutne MAT_CO
FEM3D 90 CRS 形式への変換 :MAT_CO!C!C***!C*** MAT_CO!C***!C subroutne MAT_CO use pfem_utl mplct REAL* (A-HO-Z) allocate (nde(0:)) nde= 0 do = nde()= nde(-) + ILU() PLU= nde() allocate (tem(plu)) do = do k= ILU() kk = k + nde(-) tem(kk)= IALU(k) tem に から始まる節点番号を記憶 deallocate (ILU IALU) end subroutne MAT_CO
FEM3D 9 CRS 形式への変換 :MAT_CO!C!C***!C*** MAT_CO!C***!C subroutne MAT_CO use pfem_utl mplct REAL* (A-HO-Z) allocate (nde(0:)) nde= 0 do = nde()= nde(-) + ILU() PLU= nde() allocate (tem(plu)) do = do k= ILU() kk = k + nde(-) tem(kk)= IALU(k) deallocate (ILU IALU) end subroutne MAT_CO これらはもはや不要
FEM3D 92 全体処理 program heat3d use solver use pfem_utl mplct REAL*(A-HO-Z) call IPUT_CTL call IPUT_GRID call MAT_CO0 call MAT_CO call MAT_ASS_MAI call MAT_ASS_BC call SOLVE call OUTPUT_UCD end program heat3d
FEM3D 93 MAT_ASS_MAI: 全体構成 do kpn= 2 ガウス積分点番号 ( 方向 ) do jpn= 2 ガウス積分点番号 ( 方向 ) do pn= 2 ガウス積分点番号 ( 方向 ) ガウス積分点 ( 個 ) における形状関数 およびその 自然座標系 における微分の算出 do cel= ICELTOT 要素ループ 節点の座標から ガウス積分点における 形状関数の 全体座標系 における微分 およびヤコビアンを算出 (JACOBI) do e= 局所節点番号 do je= 局所節点番号全体節点番号 :p jp A pjp のtemにおけるアドレス :kk j e do kpn= 2 ガウス積分点番号 ( 方向 ) do jpn= 2 ガウス積分点番号 ( 方向 ) do pn= 2 ガウス積分点番号 ( 方向 ) 要素積分 要素行列成分計算 全体行列への足しこみ e
FEM3D 94 係数行列 :MAT_ASS_MAI(/6)!C!C***!C*** MAT_ASS_MAI!C***!C subroutne MAT_ASS_MAI use pfem_utl mplct REAL* (A-HO-Z) nteger(knd=knt) dmenson( ) :: nodlocal allocate (AMAT(PLU)) allocate (B() D() X()) AMAT= 0.d0 係数行列 ( 非零非対角成分 ) B= 0.d0 右辺ベクトル X= 0.d0 未知数ベクトル D= 0.d0 係数行列 ( 対角成分 ) WEI()= +.0000000000D+00 WEI(2)= +.0000000000D+00 POS()= -0.5773502692D+00 POS(2)= +0.5773502692D+00
FEM3D 95 係数行列 :MAT_ASS_MAI(/6)!C!C***!C*** MAT_ASS_MAI!C***!C subroutne MAT_ASS_MAI use pfem_utl mplct REAL* (A-HO-Z) nteger(knd=knt) dmenson( ) :: nodlocal allocate (AMAT(PLU)) allocate (B() D() X()) AMAT= 0.d0 B= 0.d0 X= 0.d0 D= 0.d0 WEI()= +.0000000000D+00 WEI(2)= +.0000000000D+00 POS()= -0.5773502692D+00 POS(2)= +0.5773502692D+00 POS: WEI: 積分点座標重み係数
FEM3D 96 係数行列 :MAT_ASS_MAI(2/6)!C!C-- IIT.!C PQ - st-order dervatve of shape functon b QSI!C PE - st-order dervatve of shape functon b ETA!C PT - st-order dervatve of shape functon b ZET!C do kp= 2 do jp= 2 do p= 2 QP=.d0 + POS(p) QM=.d0 - POS(p) EP=.d0 + POS(jp) EM=.d0 - POS(jp) TP=.d0 + POS(kp) TM=.d0 - POS(kp) SHAPE(pjpkp)= Oth * QM * EM * TM SHAPE(pjpkp2)= Oth * QP * EM * TM SHAPE(pjpkp3)= Oth * QP * EP * TM SHAPE(pjpkp4)= Oth * QM * EP * TM SHAPE(pjpkp5)= Oth * QM * EM * TP SHAPE(pjpkp6)= Oth * QP * EM * TP SHAPE(pjpkp7)= Oth * QP * EP * TP
FEM3D 97 係数行列 :MAT_ASS_MAI(2/6)!C!C-- IIT.!C PQ - st-order dervatve of shape functon b QSI!C PE - st-order dervatve of shape functon b ETA!C PT - st-order dervatve of shape functon b ZET!C do kp= 2 do jp= 2 do p= 2 QP=.d0 + POS(p) QM=.d0 - POS(p) EP=.d0 + POS(jp) EM=.d0 - POS(jp) TP=.d0 + POS(kp) TM=.d0 - POS(kp) SHAPE(pjpkp)= Oth * QM * EM * TM SHAPE(pjpkp2)= Oth * QP * EM * TM SHAPE(pjpkp3)= Oth * QP * EP * TM SHAPE(pjpkp4)= Oth * QM * EP * TM SHAPE(pjpkp5)= Oth * QM * EM * TP SHAPE(pjpkp6)= Oth * QP * EM * TP SHAPE(pjpkp7)= Oth * QP * EP * TP QP EP TP k QM j EM j j TMk k k
FEM3D 9 係数行列 :MAT_ASS_MAI(2/6)!C!C-- IIT.!C PQ - st-order dervatve of shape functon b QSI!C PE - st-order dervatve of shape functon b ETA!C PT - st-order dervatve of shape functon b ZET!C do kp= 2 do jp= 2 do p= 2 QP=.d0 + POS(p) QM=.d0 - POS(p) EP=.d0 + POS(jp) EM=.d0 - POS(jp) TP=.d0 + POS(kp) TM=.d0 - POS(kp) SHAPE(pjpkp)= Oth * QM * EM * TM SHAPE(pjpkp2)= Oth * QP * EM * TM SHAPE(pjpkp3)= Oth * QP * EP * TM SHAPE(pjpkp4)= Oth * QM * EP * TM SHAPE(pjpkp5)= Oth * QM * EM * TP SHAPE(pjpkp6)= Oth * QP * EM * TP SHAPE(pjpkp7)= Oth * QP * EP * TP 7 5 6 4 2 3
FEM3D 99 係数行列 :MAT_ASS_MAI(2/6)!C!C-- IIT.!C PQ - st-order dervatve of shape functon b QSI!C PE - st-order dervatve of shape functon b ETA!C PT - st-order dervatve of shape functon b ZET!C do kp= 2 do jp= 2 do p= 2 QP=.d0 + POS(p) QM=.d0 - POS(p) EP=.d0 + POS(jp) EM=.d0 - POS(jp) TP=.d0 + POS(kp) TM=.d0 - POS(kp) SHAPE(pjpkp)= Oth * QM * EM * TM SHAPE(pjpkp2)= Oth * QP * EM * TM SHAPE(pjpkp3)= Oth * QP * EP * TM SHAPE(pjpkp4)= Oth * QM * EP * TM SHAPE(pjpkp5)= Oth * QM * EM * TP SHAPE(pjpkp6)= Oth * QP * EM * TP SHAPE(pjpkp7)= Oth * QP * EP * TP ) ( ) ( ) ( ) ( 4 3 2 ) ( ) ( ) ( ) ( 7 6 5
FEM3D 00 係数行列 :MAT_ASS_MAI(3/6) PQ(jpkp)= - Oth * EM * TM PQ(jpkp2)= + Oth * EM * TM PQ(jpkp3)= + Oth * EP * TM PQ(jpkp4)= - Oth * EP * TM PQ(jpkp5)= - Oth * EM * TP PQ(jpkp6)= + Oth * EM * TP PQ(jpkp7)= + Oth * EP * TP PQ(jpkp)= - Oth * EP * TP PE(pkp)= - Oth * QM * TM PE(pkp2)= - Oth * QP * TM PE(pkp3)= + Oth * QP * TM PE(pkp4)= + Oth * QM * TM PE(pkp5)= - Oth * QM * TP PE(pkp6)= - Oth * QP * TP PE(pkp7)= + Oth * QP * TP PE(pkp)= + Oth * QM * TP PT(pjp)= - Oth * QM * EM PT(pjp2)= - Oth * QP * EM PT(pjp3)= - Oth * QP * EP PT(pjp4)= - Oth * QM * EP PT(pjp5)= + Oth * QM * EM PT(pjp6)= + Oth * QP * EM PT(pjp7)= + Oth * QP * EP PT(pjp)= + Oth * QM * EP do cel= ICELTOT COD0= COD n= ICELOD(cel) n2= ICELOD(cel2) n3= ICELOD(cel3) n4= ICELOD(cel4) n5= ICELOD(cel5) n6= ICELOD(cel6) n7= ICELOD(cel7) n= ICELOD(cel) ( j k ) PQ( j k) l ( j k ) PE( k) l ( j k ) PT j) l ( j ) ( k ( j k ) 2 ( j k ) 3 ( j k ) 3 ( j k ) j j j j における形状関数の一階微分 k k k k
FEM3D 0 係数行列 :MAT_ASS_MAI(3/6) PQ(jpkp)= - Oth * EM * TM PQ(jpkp2)= + Oth * EM * TM PQ(jpkp3)= + Oth * EP * TM PQ(jpkp4)= - Oth * EP * TM PQ(jpkp5)= - Oth * EM * TP PQ(jpkp6)= + Oth * EM * TP PQ(jpkp7)= + Oth * EP * TP PQ(jpkp)= - Oth * EP * TP PE(pkp)= - Oth * QM * TM PE(pkp2)= - Oth * QP * TM PE(pkp3)= + Oth * QP * TM PE(pkp4)= + Oth * QM * TM PE(pkp5)= - Oth * QM * TP PE(pkp6)= - Oth * QP * TP PE(pkp7)= + Oth * QP * TP PE(pkp)= + Oth * QM * TP PT(pjp)= - Oth * QM * EM PT(pjp2)= - Oth * QP * EM PT(pjp3)= - Oth * QP * EP PT(pjp4)= - Oth * QM * EP PT(pjp5)= + Oth * QM * EM PT(pjp6)= + Oth * QP * EM PT(pjp7)= + Oth * QP * EP PT(pjp)= + Oth * QM * EP do cel= ICELTOT COD0= COD n= ICELOD(cel) n2= ICELOD(cel2) n3= ICELOD(cel3) n4= ICELOD(cel4) n5= ICELOD(cel5) n6= ICELOD(cel6) n7= ICELOD(cel7) n= ICELOD(cel) 7 5 6 4 3 2
FEM3D 02 係数行列 :MAT_ASS_MAI(4/6) nodlocal()= n nodlocal(2)= n2 nodlocal(3)= n3 nodlocal(4)= n4 nodlocal(5)= n5 nodlocal(6)= n6 nodlocal(7)= n7 nodlocal()= n 節点の節点番号 X= XYZ(n) X2= XYZ(n2) X3= XYZ(n3) X4= XYZ(n4) X5= XYZ(n5) X6= XYZ(n6) X7= XYZ(n7) X= XYZ(n) Y= XYZ(n2) Y2= XYZ(n22) Y3= XYZ(n32) Y4= XYZ(n42) Y5= XYZ(n52) Y6= XYZ(n62) Y7= XYZ(n72) Y= XYZ(n2) QVC= Oth * (X+X2+X3+X4+X5+X6+X7+X+ & Y+Y2+Y3+Y4+Y5+Y6+Y7+Y) Z= XYZ(n3) Z2= XYZ(n23) Z3= XYZ(n33) Z4= XYZ(n43) Z5= XYZ(n53) Z6= XYZ(n63) Z7= XYZ(n73) Z= XYZ(n3) 7 5 6 4 2 3 call JACOBI (DETJ PQ PE PT PX PY PZ & & X X2 X3 X4 X5 X6 X7 X & & Y Y2 Y3 Y4 Y5 Y6 Y7 Y & & Z Z2 Z3 Z4 Z5 Z6 Z7 Z )
FEM3D 03 係数行列 :MAT_ASS_MAI(4/6) nodlocal()= n nodlocal(2)= n2 nodlocal(3)= n3 nodlocal(4)= n4 nodlocal(5)= n5 nodlocal(6)= n6 nodlocal(7)= n7 nodlocal()= n X= XYZ(n) X2= XYZ(n2) X3= XYZ(n3) X4= XYZ(n4) X5= XYZ(n5) X6= XYZ(n6) X7= XYZ(n7) X= XYZ(n) Y= XYZ(n2) Y2= XYZ(n22) Y3= XYZ(n32) Y4= XYZ(n42) Y5= XYZ(n52) Y6= XYZ(n62) Y7= XYZ(n72) Y= XYZ(n2) QVC= Oth * (X+X2+X3+X4+X5+X6+X7+X+ & Y+Y2+Y3+Y4+Y5+Y6+Y7+Y) Z= XYZ(n3) Z2= XYZ(n23) Z3= XYZ(n33) Z4= XYZ(n43) Z5= XYZ(n53) Z6= XYZ(n63) Z7= XYZ(n73) Z= XYZ(n3) 節点の X 座標 節点の Y 座標 節点の Z 座標 7 5 6 4 2 3 call JACOBI (DETJ PQ PE PT PX PY PZ & & X X2 X3 X4 X5 X6 X7 X & & Y Y2 Y3 Y4 Y5 Y6 Y7 Y & & Z Z2 Z3 Z4 Z5 Z6 Z7 Z )
FEM3D 04 係数行列 :MAT_ASS_MAI(4/6) nodlocal()= n nodlocal(2)= n2 nodlocal(3)= n3 nodlocal(4)= n4 nodlocal(5)= n5 nodlocal(6)= n6 nodlocal(7)= n7 nodlocal()= n X= XYZ(n) X2= XYZ(n2) X3= XYZ(n3) X4= XYZ(n4) X5= XYZ(n5) X6= XYZ(n6) X7= XYZ(n7) X= XYZ(n) Y= XYZ(n2) Y2= XYZ(n22) Y3= XYZ(n32) Y4= XYZ(n42) Y5= XYZ(n52) Y6= XYZ(n62) Y7= XYZ(n72) Y= XYZ(n2) QVC= Oth * (X+X2+X3+X4+X5+X6+X7+X+ & Y+Y2+Y3+Y4+Y5+Y6+Y7+Y) Z= XYZ(n3) Z2= XYZ(n23) Z3= XYZ(n33) Z4= XYZ(n43) Z5= XYZ(n53) Z6= XYZ(n63) Z7= XYZ(n73) Z= XYZ(n3) 節点の X 座標 節点の Y 座標 T Q 2 5 6 T 4 z QVOL C C T z z 7 3 Q z 0 体積当たり発熱量は位置 ( メッシュの中心の座標 c c ) に依存 call JACOBI (DETJ PQ PE PT PX PY PZ & & X X2 X3 X4 X5 X6 X7 X & & Y Y2 Y3 Y4 Y5 Y6 Y7 Y & & Z Z2 Z3 Z4 Z5 Z6 Z7 Z )
FEM3D 05 係数行列 :MAT_ASS_MAI(4/6) nodlocal()= n nodlocal(2)= n2 nodlocal(3)= n3 nodlocal(4)= n4 nodlocal(5)= n5 nodlocal(6)= n6 nodlocal(7)= n7 nodlocal()= n X= XYZ(n) X2= XYZ(n2) X3= XYZ(n3) X4= XYZ(n4) X5= XYZ(n5) X6= XYZ(n6) X7= XYZ(n7) X= XYZ(n) Y= XYZ(n2) Y2= XYZ(n22) Y3= XYZ(n32) Y4= XYZ(n42) Y5= XYZ(n52) Y6= XYZ(n62) Y7= XYZ(n72) Y= XYZ(n2) QVC= Oth * (X+X2+X3+X4+X5+X6+X7+X+ & Y+Y2+Y3+Y4+Y5+Y6+Y7+Y) Z= XYZ(n3) Z2= XYZ(n23) Z3= XYZ(n33) Z4= XYZ(n43) Z5= XYZ(n53) Z6= XYZ(n63) Z7= XYZ(n73) Z= XYZ(n3) T Q 2 5 6 T 4 z QVOL C C QVC C C T z z 7 3 Q z 0 call JACOBI (DETJ PQ PE PT PX PY PZ & & X X2 X3 X4 X5 X6 X7 X & & Y Y2 Y3 Y4 Y5 Y6 Y7 Y & & Z Z2 Z3 Z4 Z5 Z6 Z7 Z )
FEM3D 06 係数行列 :MAT_ASS_MAI(4/6) nodlocal()= n nodlocal(2)= n2 nodlocal(3)= n3 nodlocal(4)= n4 nodlocal(5)= n5 nodlocal(6)= n6 nodlocal(7)= n7 nodlocal()= n X= XYZ(n) X2= XYZ(n2) X3= XYZ(n3) X4= XYZ(n4) X5= XYZ(n5) X6= XYZ(n6) X7= XYZ(n7) X= XYZ(n) Y= XYZ(n2) Y2= XYZ(n22) Y3= XYZ(n32) Y4= XYZ(n42) Y5= XYZ(n52) Y6= XYZ(n62) Y7= XYZ(n72) Y= XYZ(n2) QVC= Oth * (X+X2+X3+X4+X5+X6+X7+X+ & Y+Y2+Y3+Y4+Y5+Y6+Y7+Y) Z= XYZ(n3) Z2= XYZ(n23) Z3= XYZ(n33) Z4= XYZ(n43) Z5= XYZ(n53) Z6= XYZ(n63) Z7= XYZ(n73) Z= XYZ(n3) call JACOBI (DETJ PQ PE PT PX PY PZ & & X X2 X3 X4 X5 X6 X7 X & & Y Y2 Y3 Y4 Y5 Y6 Y7 Y & & Z Z2 Z3 Z4 Z5 Z6 Z7 Z )
FEM3D 07 JACOBI(/4) subroutne JACOBI (DETJ PQ PE PT PX PY PZ & & X X2 X3 X4 X5 X6 X7 X Y Y2 Y3 Y4 Y5 Y6 Y7 Y & & Z Z2 Z3 Z4 Z5 Z6 Z7 Z )!C!C calculates JACOBIA & IVERSE JACOBIA!C d/d d/d & d/dz!c mplct REAL* (A-HO-Z) dmenson DETJ(222) dmenson PQ(22) PE(22) PT(22) dmenson PX(222) PY(222) PZ(222) do kp= 2 do jp= 2 do p= 2 PX(pjpkp)=0.d0 PX(pjpkp2)=0.d0 PX(pjpkp3)=0.d0 PX(pjpkp4)=0.d0 PX(pjpkp5)=0.d0 PX(pjpkp6)=0.d0 PX(pjpkp7)=0.d0 PX(pjpkp)=0.d0 PY(pjpkp)=0.d0 PY(pjpkp2)=0.d0 PY(pjpkp3)=0.d0 PY(pjpkp4)=0.d0 PY(pjpkp5)=0.d0 PY(pjpkp6)=0.d0 PY(pjpkp7)=0.d0 PY(pjpkp)=0.d0 PZ(pjpkp)=0.d0 PZ(pjpkp2)=0.d0 PZ(pjpkp3)=0.d0 PZ(pjpkp4)=0.d0 PZ(pjpkp5)=0.d0 PZ(pjpkp6)=0.d0 PZ(pjpkp7)=0.d0 PZ(pjpkp)=0.d0 l l l 入力 z l ~ 出力 l l l det J z 各ガウス積分点 (pjpkp) における値 l l l
FEM3D 0 自然座標系における偏微分 (/4) 偏微分の公式より以下のようになる : z z z z z z ) ( ) ( ) ( z は定義より簡単に求められるがを実際の計算で使用する
FEM3D 09 自然座標系における偏微分 (2/4) マトリックス表示すると : z J z z z z J : ヤコビのマトリクス (Jacob matr Jacoban) 33 32 3 23 22 2 3 2 J J J J J J J J J J
FEM3D 0 自然座標系における偏微分 (3/4) の定義より簡単に求められる z z z J J J z z z J J J z z z J J J 33 32 3 23 22 2 3 2
FEM3D!C!C== DETERMIAT of the JACOBIA JACOBI(2/4) dxdq = & & + PQ(jpkp) * X + PQ(jpkp2) * X2 & & + PQ(jpkp3) * X3 + PQ(jpkp4) * X4 & & + PQ(jpkp5) * X5 + PQ(jpkp6) * X6 & & + PQ(jpkp7) * X7 + PQ(jpkp) * X dydq = & & + PQ(jpkp) * Y + PQ(jpkp2) * Y2 & & + PQ(jpkp3) * Y3 + PQ(jpkp4) * Y4 & J J2 & + PQ(jpkp5) * Y5 + PQ(jpkp6) * Y6 & & + PQ(jpkp7) * Y7 + PQ(jpkp) * Y dzdq = & J & + PQ(jpkp) * Z + PQ(jpkp2) * Z2 & J 2 J 22 & + PQ(jpkp3) * Z3 + PQ(jpkp4) * Z4 & & + PQ(jpkp5) * Z5 + PQ(jpkp6) * Z6 & J3 J32 & + PQ(jpkp7) * Z7 + PQ(jpkp) * Z dxde = & & + PE(pkp) * X + PE(pkp2) * X2 & & + PE(pkp3) * X3 + PE(pkp4) * X4 & & + PE(pkp5) * X5 + PE(pkp6) * X6 & & + PE(pkp7) * X7 + PE(pkp) * X dyde = & dxdq J & + PE(pkp) * Y + PE(pkp2) * Y2 & & + PE(pkp3) * Y3 + PE(pkp4) * Y4 & & + PE(pkp5) * Y5 + PE(pkp6) * Y6 & & + PE(pkp7) * Y7 + PE(pkp) * Y dydq J2 dzde = & & + PE(pkp) * Z + PE(pkp2) * Z2 & & + PE(pkp3) * Z3 + PE(pkp4) * Z4 & & + PE(pkp5) * Z5 + PE(pkp6) * Z6 & z & + PE(pkp7) * Z7 + PE(pkp) * Z dzdq J3 J J J 3 23 33
FEM3D 2 JACOBI(3/4) dxdt = & & + PT(pjp) * X + PT(pjp2) * X2 & & + PT(pjp3) * X3 + PT(pjp4) * X4 & & + PT(pjp5) * X5 + PT(pjp6) * X6 & & + PT(pjp7) * X7 + PT(pjp) * X dydt = & & + PT(pjp) * Y + PT(pjp2) * Y2 & & + PT(pjp3) * Y3 + PT(pjp4) * Y4 & & + PT(pjp5) * Y5 + PT(pjp6) * Y6 & & + PT(pjp7) * Y7 + PT(pjp) * Y dzdt = & & + PT(pjp) * Z + PT(pjp2) * Z2 & & + PT(pjp3) * Z3 + PT(pjp4) * Z4 & & + PT(pjp5) * Z5 + PT(pjp6) * Z6 & & + PT(pjp7) * Z7 + PT(pjp) * Z DETJ(pjpkp)= dxdq*(dyde*dzdt-dzde*dydt) + & & dydq*(dzde*dxdt-dxde*dzdt) + & & dzdq*(dxde*dydt-dyde*dxdt)!c!c== IVERSE JACOBIA coef=.d0 / DETJ(pjpkp) a= coef * ( dyde*dzdt - dzde*dydt ) a2= coef * ( dzdq*dydt - dydq*dzdt ) a3= coef * ( dydq*dzde - dzdq*dyde ) a2= coef * ( dzde*dxdt - dxde*dzdt ) a22= coef * ( dxdq*dzdt - dzdq*dxdt ) a23= coef * ( dzdq*dxde - dxdq*dzde ) J J J J 2 3 J J J 2 22 32 J J J 3 23 33 a3= coef * ( dxde*dydt - dyde*dxdt ) a32= coef * ( dydq*dxdt - dxdq*dydt ) a33= coef * ( dxdq*dyde - dydq*dxde ) DETJ(pjpkp)= dabs(detj(pjpkp))
FEM3D 3 自然座標系における偏微分 (4/4) 従って下記のように偏微分を計算できる ヤコビアン (3 3 行列 ) の逆行列を求める J z z z z
FEM3D 4 JACOBI(3/4) dxdt = & & + PT(pjp) * X + PT(pjp2) * X2 & & + PT(pjp3) * X3 + PT(pjp4) * X4 & & + PT(pjp5) * X5 + PT(pjp6) * X6 & & + PT(pjp7) * X7 + PT(pjp) * X dydt = & & + PT(pjp) * Y + PT(pjp2) * Y2 & & + PT(pjp3) * Y3 + PT(pjp4) * Y4 & & + PT(pjp5) * Y5 + PT(pjp6) * Y6 & & + PT(pjp7) * Y7 + PT(pjp) * Y dzdt = & & + PT(pjp) * Z + PT(pjp2) * Z2 & & + PT(pjp3) * Z3 + PT(pjp4) * Z4 & & + PT(pjp5) * Z5 + PT(pjp6) * Z6 & & + PT(pjp7) * Z7 + PT(pjp) * Z DETJ(pjpkp)= dxdq*(dyde*dzdt-dzde*dydt) + & & dydq*(dzde*dxdt-dxde*dzdt) + & & dzdq*(dxde*dydt-dyde*dxdt)!c!c== IVERSE JACOBIA coef=.d0 / DETJ(pjpkp) a= coef * ( dyde*dzdt - dzde*dydt ) a2= coef * ( dzdq*dydt - dydq*dzdt ) a3= coef * ( dydq*dzde - dzdq*dyde ) a2= coef * ( dzde*dxdt - dxde*dzdt ) a22= coef * ( dxdq*dzdt - dzdq*dxdt ) a23= coef * ( dzdq*dxde - dxdq*dzde ) a3= coef * ( dxde*dydt - dyde*dxdt ) a32= coef * ( dydq*dxdt - dxdq*dydt ) a33= coef * ( dxdq*dyde - dydq*dxde ) J a a a 2 3 a a a 2 22 32 a a a 3 23 33 DETJ(pjpkp)= dabs(detj(pjpkp))
FEM3D 5 JACOBI(4/4)!C!C== set the d/dx d/dy & d/dz components PX(pjpkp)= a*pq(jpkp) + a2*pe(pkp) + a3*pt(pjp) PX(pjpkp2)= a*pq(jpkp2) + a2*pe(pkp2) + a3*pt(pjp2) PX(pjpkp3)= a*pq(jpkp3) + a2*pe(pkp3) + a3*pt(pjp3) PX(pjpkp4)= a*pq(jpkp4) + a2*pe(pkp4) + a3*pt(pjp4) PX(pjpkp5)= a*pq(jpkp5) + a2*pe(pkp5) + a3*pt(pjp5) PX(pjpkp6)= a*pq(jpkp6) + a2*pe(pkp6) + a3*pt(pjp6) PX(pjpkp7)= a*pq(jpkp7) + a2*pe(pkp7) + a3*pt(pjp7) PX(pjpkp)= a*pq(jpkp) + a2*pe(pkp) + a3*pt(pjp) PY(pjpkp)= a2*pq(jpkp) + a22*pe(pkp) + a23*pt(pjp) PY(pjpkp2)= a2*pq(jpkp2) + a22*pe(pkp2) + a23*pt(pjp2) PY(pjpkp3)= a2*pq(jpkp3) + a22*pe(pkp3) + a23*pt(pjp3) PY(pjpkp4)= a2*pq(jpkp4) + a22*pe(pkp4) + a23*pt(pjp4) PY(pjpkp5)= a2*pq(jpkp5) + a22*pe(pkp5) + a23*pt(pjp5) PY(pjpkp6)= a2*pq(jpkp6) + a22*pe(pkp6) + a23*pt(pjp6) PY(pjpkp7)= a2*pq(jpkp7) + a22*pe(pkp7) + a23*pt(pjp7) PY(pjpkp)= a2*pq(jpkp) + a22*pe(pkp) + a23*pt(pjp) PZ(pjpkp)= a3*pq(jpkp) + a32*pe(pkp) + a33*pt(pjp) PZ(pjpkp2)= a3*pq(jpkp2) + a32*pe(pkp2) + a33*pt(pjp2) PZ(pjpkp3)= a3*pq(jpkp3) + a32*pe(pkp3) + a33*pt(pjp3) PZ(pjpkp4)= a3*pq(jpkp4) + a32*pe(pkp4) + a33*pt(pjp4) PZ(pjpkp5)= a3*pq(jpkp5) + a32*pe(pkp5) + a33*pt(pjp5) PZ(pjpkp6)= a3*pq(jpkp6) + a32*pe(pkp6) + a33*pt(pjp6) PZ(pjpkp7)= a3*pq(jpkp7) + a32*pe(pkp7) + a33*pt(pjp7) PZ(pjpkp)= a3*pq(jpkp) + a32*pe(pkp) + a33*pt(pjp) a a a a a a a a a z z z z 33 32 3 23 22 2 3 2
FEM3D 6 係数行列 :MAT_ASS_MAI(5/6)!C!C== COSTRUCT the GLOBAL MATRIX do e= p = nodlocal(e) do je= jp = nodlocal(je) kk= 0 f (jp.ne.p) then S= nde(p-) + E= nde(p ) do k= S E f ( tem(k).eq.jp ) then kk= k et endf endf 5 6 4 7 3 全体行列の非対角成分 A p jp kk:tem におけるアドレス p= nodlocal(e) jp= nodlocal(je) 2
FEM3D 7 要素マトリクス : 行列 j j k j 7 5 6 4 3 2
FEM3D 係数行列 :MAT_ASS_MAI(5/6)!C!C== COSTRUCT the GLOBAL MATRIX do e= p = nodlocal(e) do je= jp = nodlocal(je) kk= 0 f (jp.ne.p) then S= nde(p-) + E= nde(p ) do k= S E f ( tem(k).eq.jp ) then kk= k et endf endf 要素マトリクス ( e ~j e ) 全体マトリクス ( p ~j p ) の関係 kk:tem におけるアドレス 7 5 6 4 3 2
FEM3D 9 係数行列 :MAT_ASS_MAI(6/6) QV0 = 0.d0 COEFj= 0.d0 do kpn= 2 do jpn= 2 do pn= 2 coef= dabs(detj(pnjpnkpn))*wei(pn)*wei(jpn)*wei(kpn) PX= PX(pnjpnkpne) PY= PY(pnjpnkpne) PZ= PZ(pnjpnkpne) PXj= PX(pnjpnkpnje) PYj= PY(pnjpnkpnje) PZj= PZ(pnjpnkpnje) COEFj= COEFj + coef * COD0 * & (PX*PXj+PY*PYj+PZ*PZj) SH= SHAPE(pnjpnkpne) QV0= QV0 + SH * QVOL * coef f (jp.eq.p) then D(p)= D(p) + COEFj B(p)= B(p) + QV0*QVC else AMAT(kk)= AMAT(kk) + COEFj endf return end j j z z j det J ddd
FEM3D 20 係数行列 :MAT_ASS_MAI(6/6) QV0 = 0.d0 COEFj= 0.d0 do kpn= 2 do jpn= 2 do pn= 2 coef= dabs(detj(pnjpnkpn))*wei(pn)*wei(jpn)*wei(kpn) PX= PX(pnjpnkpne) PY= PY(pnjpnkpne) PZ= PZ(pnjpnkpne) PXj= PX(pnjpnkpnje) PYj= PY(pnjpnkpnje) PZj= PZ(pnjpnkpnje) COEFj= COEFj + coef * COD0 * & (PX*PXj+PY*PYj+PZ*PZj) return end SH= SHAPE(pnjpnkpne) QV0= QV0 + SH * QVOL * coef f (jp.eq.p) then D(p)= D(p) + COEFj B(p)= B(p) + QV0*QVC else AMAT(kk)= AMAT(kk) + COEFj endf j I j L M j z f ( ) ddd W W j Wk f ( j k ) k z j det J ddd
FEM3D 2 係数行列 :MAT_ASS_MAI(6/6) QV0 = 0.d0 COEFj= 0.d0 do kpn= 2 do jpn= 2 do pn= 2 coef= dabs(detj(pnjpnkpn))*wei(pn)*wei(jpn)*wei(kpn) PX= PX(pnjpnkpne) PY= PY(pnjpnkpne) PZ= PZ(pnjpnkpne) coef W W W det j k J j k PXj= PX(pnjpnkpnje) PYj= PY(pnjpnkpnje) PZj= PZ(pnjpnkpnje) COEFj= COEFj + coef * COD0 * & (PX*PXj+PY*PYj+PZ*PZj) return end SH= SHAPE(pnjpnkpne) QV0= QV0 + SH * QVOL * coef f (jp.eq.p) then D(p)= D(p) + COEFj B(p)= B(p) + QV0*QVC else AMAT(kk)= AMAT(kk) + COEFj endf j I j L M j z f ( ) ddd W W j Wk f ( j k ) k z j det J ddd
FEM3D 22 係数行列 :MAT_ASS_MAI(6/6) QV0 = 0.d0 COEFj= 0.d0 do kpn= 2 do jpn= 2 do pn= 2 coef= dabs(detj(pnjpnkpn))*wei(pn)*wei(jpn)*wei(kpn) PX= PX(pnjpnkpne) PY= PY(pnjpnkpne) PZ= PZ(pnjpnkpne) PXj= PX(pnjpnkpnje) PYj= PY(pnjpnkpnje) PZj= PZ(pnjpnkpnje) COEFj= COEFj + coef * COD0 * & (PX*PXj+PY*PYj+PZ*PZj) SH= SHAPE(pnjpnkpne) QV0= QV0 + SH * QVOL * coef f (jp.eq.p) then D(p)= D(p) + COEFj B(p)= B(p) + QV0*QVC else AMAT(kk)= AMAT(kk) + COEFj endf return end j k j j
FEM3D 23 係数行列 :MAT_ASS_MAI(6/6) QV0 = 0.d0 COEFj= 0.d0 do kpn= 2 do jpn= 2 do pn= 2 coef= dabs(detj(pnjpnkpn))*wei(pn)*wei(jpn)*wei(kpn) PX= PX(pnjpnkpne) PY= PY(pnjpnkpne) PZ= PZ(pnjpnkpne) PXj= PX(pnjpnkpnje) PYj= PY(pnjpnkpnje) PZj= PZ(pnjpnkpnje) COEFj= COEFj + coef * COD0 * & (PX*PXj+PY*PYj+PZ*PZj) SH= SHAPE(pnjpnkpne) QV0= QV0 + SH * QVOL * coef f (jp.eq.p) then D(p)= D(p) + COEFj B(p)= B(p) + QV0*QVC else AMAT(kk)= AMAT(kk) + COEFj endf return end k ( e) ( e) f ( e) ( e) T f Q dv Q V z QVOL C C QVC C C QV 0 QVOL f ( e) V QV 0QVC T dv
FEM3D 24 MAT_ASS_BC: 全体構成 do = 節点ループ ( ディリクレ ) 境界条件を設定する節点をマーク (IWKX) do = 節点ループ f (IWKX().eq.) then マークされた節点だったら対応する右辺ベクトル (B) の成分 対角成分 (D) の成分の修正 ( 行 列 ) do k= nde(-)+ nde() 対応する非零非対角成分 (AMAT) の成分の修正 ( 行 ) endf Z T=0@Z=z ma do = 節点ループ do k= nde (-)+ nde () f (IWKX(tem (k)).eq.) then 対応する非零非対角成分の節点がマークされていたら対応する右辺ベクトル 非零非対角成分 (AMAT) の成分の修正 ( 列 ) endf X Y X Z Y
FEM3D 25 境界条件 :MAT_ASS_BC(/2)!C!C== Z=Zma subroutne MAT_ASS_BC use pfem_utl mplct REAL* (A-HO-Z) allocate (IWKX(2)) IWKX= 0 do n= IWKX(n)= 0 b0= - do b0= ODGRPtot f (ODGRP_AME(b0).eq.'Zma') et do b= ODGRP_IDEX(b0-)+ ODGRP_IDEX(b0) n= ODGRP_ITEM(b) IWKX(n)= 節点グループ名が Zma である節点 n において : IWKX(n)= とする
FEM3D 26 境界条件 :MAT_ASS_BC(2/2) do n= f (IWKX(n).eq.) then B(n)= 0.d0 D(n)=.d0 S= nde(n-) + E= nde(n ) do k= S E AMAT(k)= 0.d0 endf!c== do n= S= nde(n-) + E= nde(n ) do k= S E f (IWKX(tem(k)).eq.) then AMAT(k)= 0.d0 endf return end
FEM3D 27 対象とする問題 : 一次元熱伝導問題 体積当たり一様発熱 Q T Q 0 =0 ( mn ) = ma 一様な : 断面積 A 熱伝導率 体積当たり一様発熱 ( 時間当たり ) QL -3 T - 境界条件 =0 := 0 ( 固定 ) T = ma : 0( 断熱 ) 復習 : 一次元問題 Q
FEM3D 2 =0 で成立する方程式 T =0 体積当たり一様発熱 Q T Q 0 =0 ( mn ) = ma 一様な : 断面積 A 熱伝導率 体積当たり一様発熱 ( 時間当たり ) QL -3 T - 境界条件 =0 :=0 ( 固定 ) T = ma : 0( 断熱 ) Q 復習 : 一次元問題
FEM3D 29 プログラム :d.f(6/6) 第一種境界条件 @=0!C!C +---------------------+!C BOUDARY CODITIOS!C +---------------------+!C===!C!C-- X=Xmn = js= IDEX(-) T =0 対角成分 = 右辺 =0 非対角成分 =0 AMAT(jS+)= 0.d0 DIAG()=.d0 RHS ()= 0.d0!C=== do k= PLU f (ITEM(k).eq.) AMAT(k)= 0.d0 復習 : 一次元問題
FEM3D 30 プログラム :d.f(6/6) 第一種境界条件 @=0!C!C +---------------------+!C BOUDARY CODITIOS!C +---------------------+!C===!C!C-- X=Xmn = js= IDEX(-) AMAT(jS+)= 0.d0 DIAG()=.d0 RHS ()= 0.d0 T =0 対角成分 = 右辺 =0 非対角成分 =0 ゼロクリア!C=== do k= PLU f (ITEM(k).eq.) AMAT(k)= 0.d0 復習 : 一次元問題
FEM3D 3 プログラム :d.f(6/6) 第一種境界条件 @=0!C!C +---------------------+!C BOUDARY CODITIOS!C +---------------------+!C===!C!C-- X=Xmn = js= IDEX(-) T =0 対角成分 = 右辺 =0 非対角成分 =0 消去 ゼロクリア AMAT(jS+)= 0.d0 DIAG()=.d0 RHS ()= 0.d0!C=== do k= PLU f (ITEM(k).eq.) AMAT(k)= 0.d0 行列の対称性を保つため 第一種境界条件を適用している節点に対応する 列 を 右辺に移項して消去する ( 今の場合は非対角成分を 0 にするだけで良い ) 復習 : 一次元問題
FEM3D 32 第一種境界条件が T 0 の場合!C!C +---------------------+!C BOUDARY CODITIOS!C +---------------------+!C=== 行列の対称性を保つため 第一種境界条件を適用している節点に対応する 列 を 右辺に移項して消去する!C!C-- X=Xmn = js= IDEX(-) Dag j j Inde[ j] Amat k Inde[ j] k Item[ k ] Rhs j AMAT(jS+)= 0.d0 DIAG()=.d0 RHS ()= PHImn!C=== do = do k= IDEX(-)+ IDEX() f (ITEM(k).eq.) then RHS ()= RHS() AMAT(k)*PHImn AMAT(k)= 0.d0 endf 復習 : 一次元問題
FEM3D 33 第一種境界条件が T 0 の場合!C!C +---------------------+!C BOUDARY CODITIOS!C +---------------------+!C=== 行列の対称性を保つため 第一種境界条件を適用している節点に対応する 列 を 右辺に移項して消去する!C!C-- X=Xmn = js= IDEX(-)!C=== Dag AMAT(jS+)= 0.d0 j DIAG()=.d0 RHS ()= PHImn Rhs j do = do k= IDEX(-)+ IDEX() f (ITEM(k).eq.) then RHS ()= RHS() AMAT(k)*PHImn AMAT(k)= 0.d0 endf j Rhs j Inde[ j] k Inde[ j] k k Amat Amat k k s s T Item[ k mn Amat s ] s k Item[ k ] where Item( k s ) 復習 : 一次元問題
FEM3D 34 境界条件 :MAT_ASS_BC(2/2) do n= f (IWKX(n).eq.) then B(n)= 0.d0 D(n)=.d0 S= nde(n-) + E= nde(n ) do k= S E AMAT(k)= 0.d0 endf IWKX(n)= となる節点に対して対角成分 = 右辺 =0 非零対角成分 =0 ゼロクリア!C== do n= S= nde(n-) + E= nde(n ) do k= S E f (IWKX(tem(k)).eq.) then AMAT(k)= 0.d0 endf return end ここでやっていることも一次元の時と全く変わらない
FEM3D 35 境界条件 :MAT_ASS_BC(2/2)!C== do n= f (IWKX(n).eq.) then B(n)= 0.d0 D(n)=.d0 S= nde(n-) + E= nde(n ) do k= S E AMAT(k)= 0.d0 endf do n= S= nde(n-) + E= nde(n ) do k= S E f (IWKX(tem(k)).eq.) then AMAT(k)= 0.d0 endf return end IWKX(n)= となる節点を非零非対角成分として有している節点に対して 右辺へ移項 当該非零非対角成分 =0 消去 ゼロクリア ここでやっていることも一次元の時と全く変わらない
FEM3D 36 test メインプログラム nput_cntl 制御データ入力 nput_grd メッシュファイル入力 fnd_node 節点探索 mat_con0 行列コネクティビティ生成 msort ソート mat_con 行列コネクティビティ生成 mat_ass_man 係数行列生成 jacob ヤコビアン計算 mat_ass_bc 境界条件処理 三次元熱伝導解 析コード heat3d の構成 solve 線形ソルバー制御 output_ucd 可視化処理 cg CG 法計算
FEM3D 37 全体処理 program heat3d use solver use pfem_utl mplct REAL*(A-HO-Z) call IPUT_CTL call IPUT_GRID call MAT_CO0 call MAT_CO call MAT_ASS_MAI call MAT_ASS_BC call SOLVE call OUTPUT_UCD end program heat3d
FEM3D 3 SOLVE module SOLVER contans subroutne SOLVE use pfem_utl use solver_cg mplct REAL* (A-HO-Z) nteger :: ERROR ICFLAG character(len=char_length) :: BUF data ICFLAG/0/!C!C +------------+!C PARAMETERs!C +------------+!C=== ITER = pfemiarra() CG 法の最大反復回数 RESID = pfemrarra() CG 法の収束判定値!C===!C!C +------------------+!C ITERATIVE solver!c +------------------+!C=== call CG & & ( PLU D AMAT nde tem B X RESID ITER ERROR )!C=== ITERactual= ITER end subroutne SOLVE end module SOLVER
FEM3D 39 前処理付き共役勾配法 Precondtoned Conjugate Gradent Method (CG) Compute r (0) = b-[a] (0) for = 2 solve [M]z (-) = r (-) - = r (-) z (-) f = p () = z (0) else 前処理 : 対角スケーリング end - = - / -2 p () = z (-) + - endf q () = [A]p () = - /p () q () () = (-) + p () r () = r (-) - q () check convergence r p (-)
FEM3D 40 対角スケーリング 点ヤコビ前処理 前処理行列として もとの行列の対角成分のみを取り出した行列を前処理行列 [M] とする 対角スケーリング 点ヤコビ (pont-jacob) 前処理 M D 0... 0 0 0 D 2 0 0......... D 0 0 0 0 0... 0 D solve [M]z (-) = r (-) という場合に逆行列を簡単に求めることができる
FEM3D 4 module solver_cg contans CG ソルバー (/6) subroutne CG & & ( PLU D AMAT nde tem B X RESID ITER ERROR) mplct REAL*(A-HO-Z) nclude 'precson.nc' nteger(knd=knt ) ntent(n):: PLU nteger(knd=knt ) ntent(nout):: ITER ERROR real (knd=kreal) ntent(nout):: RESID real(knd=kreal) dmenson() ntent(nout):: B X D real(knd=kreal) dmenson(plu) ntent(nout):: AMAT nteger(knd=knt ) dmenson(0: )ntent(n) :: nde nteger(knd=knt ) dmenson(plu)ntent(n) :: tem real(knd=kreal) dmenson(::) allocatable :: WW nteger(knd=knt) parameter :: R= nteger(knd=knt) parameter :: Z= 2 nteger(knd=knt) parameter :: Q= 2 nteger(knd=knt) parameter :: P= 3 nteger(knd=knt) parameter :: DD= 4 nteger(knd=knt ) :: MAXIT real (knd=kreal) :: TOL W SS
FEM3D 42 module solver_cg contans CG ソルバー (/6) WW()= WW(R) {r} WW(2)= WW(Z) {z} WW(2)= mplct REAL*(A-HO-Z) WW(Q) {q} nclude 'precson.nc' WW(3)= WW(P) {p} WW(4)= WW(DD) /{D} subroutne CG & & ( PLU D AMAT nde tem B X RESID ITER ERROR) nteger(knd=knt ) ntent(n):: PLU nteger(knd=knt ) ntent(nout):: ITER ERROR real (knd=kreal) ntent(nout):: RESID real(knd=kreal) dmenson() ntent(nout):: B X D real(knd=kreal) dmenson(plu) ntent(nout):: AMAT nteger(knd=knt ) dmenson(0: )ntent(n) :: nde nteger(knd=knt ) dmenson(plu)ntent(n) :: tem real(knd=kreal) dmenson(::) allocatable :: WW nteger(knd=knt) parameter :: R= nteger(knd=knt) parameter :: Z= 2 nteger(knd=knt) parameter :: Q= 2 nteger(knd=knt) parameter :: P= 3 nteger(knd=knt) parameter :: DD= 4 nteger(knd=knt ) :: MAXIT real (knd=kreal) :: TOL W SS Compute r (0) = b-[a] (0) for = 2 solve [M]z (-) = r (-) - = r (-) z (-) f = p () = z (0) else - = - / -2 end p () = z (-) + - endf q () = [A]p () = - /p () q () () = (-) + p () r () = r (-) - q () check convergence r p (-)
FEM3D 43 CG ソルバー (2/6)!C!C +-------+!C IIT.!C +-------+!C=== ERROR= 0 allocate (WW(4)) MAXIT = ITER TOL = RESID X = 0.d0!C===!C +-----------------------+!C {r0}= {b} - [A]{n}!C +-----------------------+!C=== do j= WW(jDD)=.d0/D(j) WVAL= B(j) - D(j)*X(j) do k= nde(j-)+ nde(j) = tem(k) WVAL= WVAL - AMAT(k)*X() WW(jR)= WVAL WW()= WW(R) {r} WW(2)= WW(Z) {z} WW(2)= WW(Q) {q} WW(3)= WW(P) {p} WW(4)= WW(DD) /{D} 対角成分の逆数 ( 前処理用 ) その都度 除算をすると効率が悪いため 予め配列に格納
FEM3D 44 CG ソルバー (2/6)!C!C +-------+!C IIT.!C +-------+!C=== ERROR= 0 allocate (WW(4)) MAXIT = ITER TOL = RESID X = 0.d0!C===!C +-----------------------+!C {r0}= {b} - [A]{n}!C +-----------------------+!C=== do j= WW(jDD)=.d0/D(j) WVAL= B(j) - D(j)*X(j) do k= nde(j-)+ nde(j) = tem(k) WVAL= WVAL - AMAT(k)*X() WW(jR)= WVAL Compute r (0) = b-[a] (0) for = 2 solve [M]z (-) = r (-) - = r (-) z (-) f = p () = z (0) else - = - / -2 p () = z (-) + - endf q () = [A]p () = - /p () q () () = (-) + p () r () = r (-) - q () check convergence r end p (-)
FEM3D 45 CG ソルバー (3/6) BRM20= 0.d0 do = BRM20= BRM20 + B()**2 BRM2= BRM20 BRM2= b 2 あとで収束判定に使用!C=== f (BRM2.eq.0.d0) BRM2=.d0 ITER = 0 do ter= MAXIT!C!C************************************************* Conjugate Gradent Iteraton!C!C +----------------+!C {z}= [Mnv]{r}!C +----------------+!C=== do = WW(Z)= WW(R) * WW(DD)!C===
FEM3D 46 CG ソルバー (3/6)!C=== BRM20= 0.d0 do = BRM20= BRM20 + B()**2 BRM2= BRM20 f (BRM2.eq.0.d0) BRM2=.d0 ITER = 0 Compute r (0) = b-[a] (0) for = 2 solve [M]z (-) = r (-) - = r (-) z (-) f = p () = z (0) else!c!c +----------------+!C {z}= [Mnv]{r}!C +----------------+!C=== do = WW(Z)= WW(R) * WW(DD)!C=== - = - / -2 do ter= MAXIT p () = z (-) +!C - p (-)!C************************************************* Conjugate endf Gradent Iteraton end q () = [A]p () = - /p () q () () = (-) + p () r () = r (-) - q () check convergence r
FEM3D 47 CG ソルバー (4/6)!C!C +---------------+!C {RHO}= {r}{z}!c +---------------+!C=== RHO0= 0.d0 do = RHO0= RHO0 + WW(R)*WW(Z) RHO= RHO0!C===!C +-----------------------------+!C {p} = {z} f ITER=!C BETA= RHO / RHO otherwse!c +-----------------------------+!C=== f ( ITER.eq. ) then do = WW(P)= WW(Z) else BETA= RHO / RHO do = WW(P)= WW(Z) + BETA*WW(P) endf!c=== Compute r (0) = b-[a] (0) for = 2 solve [M]z (-) = r (-) - = r (-) z (-) f = p () = z (0) else - = - / -2 p () = z (-) + - endf q () = [A]p () = - /p () q () () = (-) + p () r () = r (-) - q () check convergence r end p (-)
FEM3D 4 CG ソルバー (5/6)!C +-------------+!C {q}= [A]{p}!C +-------------+!C=== do j= WVAL= D(j)*WW(jP) do k= nde(j-)+ nde(j) = tem(k) WVAL= WVAL + AMAT(k)*WW(P) WW(jQ)= WVAL!C===!C +---------------------+!C ALPHA= RHO / {p}{q}!c +---------------------+!C=== C0= 0.d0 do = C0= C0 + WW(P)*WW(Q) C= C ALPHA= RHO / C!C=== Compute r (0) = b-[a] (0) for = 2 solve [M]z (-) = r (-) - = r (-) z (-) f = p () = z (0) else - = - / -2 p () = z (-) + - endf q () = [A]p () = - /p () q () () = (-) + p () r () = r (-) - q () check convergence r end p (-)
FEM3D 49 CG ソルバー (6/6)!C!C +----------------------+!C {}= {} + ALPHA*{p}!C {r}= {r} - ALPHA*{q}!C +----------------------+!C=== do = X() = X () + ALPHA * WW(P) WW(R)= WW(R) - ALPHA * WW(Q)!C=== DRM20= 0.d0 do = DRM20= DRM20 + WW(R)**2 DRM2= DRM20 RESID= dsqrt(drm2/brm2)!c=== f ( RESID.le.TOL ) et f ( ITER.eq.MAXIT ) ERROR= -300 RHO = RHO Compute r (0) = b-[a] (0) for = 2 solve [M]z (-) = r (-) - = r (-) z (-) f = p () = z (0) else - = - / -2 p () = z (-) + - endf q () = [A]p () = - /p () q () () = (-) + p () r () = r (-) - q () check convergence r end p (-)
FEM3D 50!C!C +----------------------+!C {}= {} + ALPHA*{p}!C {r}= {r} - ALPHA*{q}!C +----------------------+!C=== do = X() = X () + ALPHA * WW(P) WW(R)= WW(R) - ALPHA * WW(Q)!C=== DRM20= 0.d0 do = DRM20= DRM20 + WW(R)**2 DRM2= DRM20 RESID= dsqrt(drm2/brm2)!c=== CG ソルバー (6/6) f ( RESID.le.TOL ) et f ( ITER.eq.MAXIT ) ERROR= -300 RHO = RHO Resd Compute r (0) = b-[a] (0) for = 2 solve [M]z (-) = r (-) - = r (-) z (-) f = p () = z (0) else - = - / -2 p () = z (-) + - endf q () = [A]p () = - /p () q () () = (-) + p () r () = r (-) - q () check convergence r end Dorm2 Borm2 r b A b b p (-) Tol