有限要素法による 三次元定常熱伝導解析プログラム C 言語編 202 年夏季集中講義中島研吾 並列計算プログラミング (66-2057) 先端計算機演習 (66-4009)
FEM3D 2 対象とする問題 : 三次元定常熱伝導 Z T T=0@Z= ma Z T 定常熱伝導 + 発熱 一様な熱伝導率 直方体 T 一辺長さの立方体 ( 六面体 ) 要素 各方向にX Y Z 個 境界条件 Q 0 X Y X Y T=0@Z= ma 体積当たり発熱量は位置 ( メッシュの中心の座標 c c ) に依存 Q QVOL C C
FEM3D 3 対象とする問題 : 三次元定常熱伝導 T T T Q 0 原点から遠い部分が高温 体積当たり発熱量は位置 ( メッシュの中心の座標 ) に依存 Q 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 4 3 - + + - 2 FEM3D 0 T T 4 4 4 ) ( ) ( ) (
FEM3D アイソパラメトリック要素 (3/3) + - + - 高次の補間関数を使えば 曲線 曲面も扱うことが可能となる そういう意味で 自然座標系 と呼んでいる Sub-Parametrc Super-Parametrc
2D 自然座標系の形状関数 (/2) 自然座標系における正方形上の内挿多項式は下式で与えられる : 各節点での条件より : 4 - + + - 3 2 FEM3D 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
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 T T T T V 0 2 2 2 2 2 2 Q T T T FEM3D 7
ガラーキン法の適用 (2/3) 三次元のグリーンの定理 前式の 2 階微分の部分に適用 ( 表面積分省略 ): これに以下を代入する : { { { { T T T T V V S dv B A B A B A ds n B A dv B B B A 2 2 2 2 2 2 dv T T T dv T T T V T T T T V FEM3D
FEM3D 9 ガラーキン法の適用 (3/3) 体積あたり発熱量の項を加えて次式が得られる : Q T T T dv Q dv 0 V この式を弱形式 (weak form) と呼ぶ 元の微分方程式では 2 階の微分が含まれていたが 上式では グリーンの定理によって 階微分に低減されている 弱形式によって近似関数 ( 形状関数 内挿関数 ) に対する要求が弱くなっている : すなわち線形関数で2 階微分の効果を記述できる 項が増えただけで 一次元と同じ V
境界条件を考慮した弱形式 : 各要素 dv dv dv k V T 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 dv V dv k j V j j j 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) 偏微分の公式より以下のようになる : ) ( ) ( ) ( は定義より簡単に求められるがを実際の計算で使用する FEM3D 24
自然座標系における偏微分 (2/4) マトリックス表示すると : J J : ヤコビのマトリクス (Jacob matr Jacoban) FEM3D 25
自然座標系における偏微分 (3/4) の定義より簡単に求められる J J J J J J J J J 33 32 3 23 22 2 3 2 FEM3D 26
自然座標系における偏微分 (4/4) 従って下記のように偏微分を計算できる ヤコビアン (3 3 行列 ) の逆行列を求める J FEM3D 27
要素単位での積分 V j j j V j j j j dv dv k 2 3 4 5 6 7 j k j FEM3D 2
自然座標系での積分 d d d J ddd dv 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 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 enddo enddo
ヒント (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= ma Z T 定常熱伝導 + 発熱 一様な熱伝導率 直方体 T 一辺長さの立方体 ( 六面体 ) 要素 各方向にX Y Z 個 境界条件 Q 0 X Y X Y T=0@Z= ma 体積当たり発熱量は位置 ( メッシュの中心の座標 c c ) に依存 Q 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= 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 QVOL C C T Q 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) FORTRA です すみません 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)= enddo enddo cou= 0 b = 2 do k= ZP j= do = XP cou= cou + = (k-)*ibodtot + (j-)*xp + IW(coub)= enddo enddo cou= 0 b = 3 k= do j= YP do = XP cou= cou + = (k-)*ibodtot + (j-)*xp + IW(coub)= enddo enddo cou= 0 b = 4 k= ZP do j= YP do = XP cou= cou + = (k-)*ibodtot + (j-)*xp + IW(coub)= enddo enddo X=Xmn の節点を IW(b) に格納 (b=yp*zp) X Z T=0@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)= enddo enddo!c=== cou= 0 b = 2 do k= ZP j= do = XP cou= cou + = (k-)*ibodtot + (j-)*xp + IW(coub)= enddo enddo cou= 0 b = 3 k= do j= YP do = XP cou= cou + = (k-)*ibodtot + (j-)*xp + IW(coub)= enddo enddo cou= 0 b = 4 k= ZP do j= YP do = XP cou= cou + = (k-)*ibodtot + (j-)*xp + IW(coub)= enddo enddo Y=Ymn の節点を IW(b2) に格納 (b=xp*zp) X Z T=0@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)= enddo enddo cou= 0 b = 2 do k= ZP j= do = XP cou= cou + = (k-)*ibodtot + (j-)*xp + IW(coub)= enddo enddo cou= 0 b = 3 k= do j= YP do = XP cou= cou + = (k-)*ibodtot + (j-)*xp + IW(coub)= enddo enddo cou= 0 b = 4 k= ZP do j= YP do = XP cou= cou + = (k-)*ibodtot + (j-)*xp + IW(coub)= enddo enddo X Z T=0@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 enddo enddo enddo 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 enddo enddo enddo 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 **/ #nclude <stdo.h> #nclude <stdlb.h> FILE* fp_log; #defne GLOBAL_VALUE_DEFIE #nclude "pfem_utl.h" //#nclude "solver.h" etern vod IPUT_CTL(); etern vod IPUT_GRID(); etern vod MAT_CO0(); etern vod MAT_CO(); etern vod MAT_ASS_MAI(); etern vod MAT_ASS_BC(); etern vod SOLVE(); etern vod OUTPUT_UCD(); nt man() { IPUT_CTL(); IPUT_GRID(); MAT_CO0(); MAT_CO(); MAT_ASS_MAI(); MAT_ASS_BC() ; SOLVE(); OUTPUT_UCD() ;
FEM3D 65 Global 変数表 :pfem_utl.h(/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 [ODGRPtot+] I 各節点グループに含まれる節点数 ( 累積 ) ODGRP_ITEM I [ODGRP_IDEX[ODG RPTOT+]] I 節点グループに含まれる節点 ODGRP_AME C0 [ODGRP_IDEX[ODG RPTOT+]] I 節点グループ名 LU I O 各節点非対角成分数 PLU I O 非対角成分総数 D R [] O 全体行列 : 対角ブロック BX R [] O 右辺ベクトル 未知数ベクトル
FEM3D 66 Global 変数表 :pfem_utl.h(2/3) 変数名種別サイズ I/O 内容 AMAT R [] O 全体行列 : 非零非対角成分 ndelu I [+] O 全体行列 : 非零非対角成分数 temlu 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.h(3/3) 変数名種別サイズ I/O 内容 Oth R I =0.25 PQ PE PT R [2][2][] O 各ガウス積分点における POS WEI R [2] O 各ガウス積分点の座標 重み係数 COL COL2 I [00] O ソート用ワーク配列 SHAPE R [2][2][2][] O 各ガウス積分点における形状関数 (=~) PX PY PZ R [2][2][2][] O 各ガウス積分点における ~ ~ DETJ R [2][2][2] O 各ガウス積分点におけるヤコビアン行列式 COD QVOL R I 熱伝導率 体積当たり発熱量係数
FEM3D 6 制御ファイル入力 :IPUT_CTL /** ** IPUT_CTL **/ #nclude <stdo.h> #nclude <stdlb.h> #nclude "pfem_utl.h" /** **/ vod IPUT_CTL() { FILE *fp; f( (fp=fopen("iput.dat""r")) == ULL){ fprntf(stdout"nput fle cannot be opened! n"); et(); fscanf(fp"%s"fname); fscanf(fp"%d"&iter); fscanf(fp "%lf %lf" &COD &QVOL); fscanf(fp "%lf" &RESID); fclose(fp); pfemrarra[0]= RESID; pfemiarra[0]= ITER; IPUT.DAT cube.0 fname 2000 ITER.0.0 COD QVOL.0e-0 RESID
FEM3D 69 メッシュ入力 :IPUT_GRID(/2) #nclude <stdo.h> #nclude <stdlb.h> #nclude "pfem_utl.h" #nclude "allocate.h" vod IPUT_GRID() { FILE *fp; nt jkkknncelse; nt TYPEIMAT; /** **/ f( (fp=fopen(fname"r")) == ULL){ fprntf(stdout"nput fle cannot be opened! n"); et(); ODE fscanf(fp"%d"&); P=; XYZ=(KREAL**)allocate_matr(seof(KREAL)3); for(=0;<;++){ for(j=0;j<3;j++){ XYZ[][j]=0.0; for(=0;<;++){ fscanf(fp"%d %lf %lf %lf"&&xyz[][0]&xyz[][]&xyz[][2]);
FEM3D 70 allocate deallocate 関数 #nclude <stdo.h> #nclude <stdlb.h> vod* allocate_vector(nt sent m) { vod *a; f ( ( a=(vod * )malloc( m * se ) ) == 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 sent mnt n) { vod **aa; nt ; f ( ( aa=(vod ** )malloc( m * seof(vod*) ) ) == ULL ) { fprntf(stdout"error:memor does not enough! aa n matr n"); et(); f ( ( aa[0]=(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) { free( aa ); allocate を FORTRA 並みに簡単にやるための関数
FEM3D 7 メッシュ入力 :IPUT_GRID(2/2) /** **/ /** **/ ELEMET fscanf(fp"%d"&iceltot); ICELOD=(KIT**)allocate_matr(seof(KIT)ICELTOT); for(=0;<iceltot;++) fscanf(fp"%d"&type); for(cel=0;cel<iceltot;cel++){ fscanf(fp"%d %d %d %d %d %d %d %d %d %d"&&imat &ICELOD[cel][0]&ICELOD[cel][]&ICELOD[cel][2]&ICELOD[cel][3] &ICELOD[cel][4]&ICELOD[cel][5]&ICELOD[cel][6]&ICELOD[cel][7]); ODE grp. nfo. fscanf(fp"%d"&odgrptot); ICELOD[][j] の中身としては から始まる通し節点番号がそのまま読み込まれている 要素番号は 0 から番号付け ODGRP_IDEX=(KIT* )allocate_vector(seof(kit)odgrptot+); ODGRP_AME =(CHAR0*)allocate_vector(seof(CHAR0)ODGRPtot); for(=0;<odgrptot+;++) ODGRP_IDEX[]=0; for(=0;<odgrptot;++) fscanf(fp"%d"&odgrp_idex[+]); nn=odgrp_idex[odgrptot]; ODGRP_ITEM=(KIT*)allocate_vector(seof(KIT)nn); for(k=0;k<odgrptot;k++){ S= ODGRP_IDEX[k]; E= ODGRP_IDEX[k+]; fscanf(fp"%s"odgrp_ame[k].name); nn= E - S; f( nn!= 0 ){ for(kk=s;kk<e;kk++) fscanf(fp"%d"&odgrp_item[kk]); 節点グループの中身も から始まる通し節点番号がそのまま読み込まれている fclose(fp);
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(/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 [ODGRPtot+] I 各節点グループに含まれる節点数 ( 累積 ) ODGRP_ITEM I [ODGRP_IDEX[ODG RPTOT+]] I 節点グループに含まれる節点 ODGRP_AME C0 [ODGRP_IDEX[ODG RPTOT+]] I 節点グループ名 LU I O 各節点非対角成分数 PLU I O 非対角成分総数 D R [] O 全体行列 : 対角ブロック BX R [] O 右辺ベクトル 未知数ベクトル
FEM3D 74 Global 変数表 :pfem_utl.h(2/3) 変数名種別サイズ I/O 内容 AMAT R [] O 全体行列 : 非零非対角成分 ndelu I [+] O 全体行列 : 非零非対角成分数 temlu 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(3/3) 変数名種別サイズ I/O 内容 Oth R I =0.25 PQ PE PT R [2][2][] O 各ガウス積分点における POS WEI R [2] O 各ガウス積分点の座標 重み係数 COL COL2 I [00] O ソート用ワーク配列 SHAPE R [2][2][2][] O 各ガウス積分点における形状関数 (=~) PX PY PZ R [2][2][2][] O 各ガウス積分点における ~ ~ DETJ R [2][2][2] 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 **/ #nclude <stdo.h> #nclude <stdlb.h> FILE* fp_log; #defne GLOBAL_VALUE_DEFIE #nclude "pfem_utl.h" //#nclude "solver.h" etern vod IPUT_CTL(); etern vod IPUT_GRID(); etern vod MAT_CO0(); etern vod MAT_CO(); etern vod MAT_ASS_MAI(); etern vod MAT_ASS_BC(); etern vod SOLVE(); etern vod OUTPUT_UCD(); nt man() { 3 4 5 6 7 9 9 0 2 4 5 6 5 6 7 2 3 2 3 4 IPUT_CTL(); IPUT_GRID(); MAT_CO0(); MAT_CO(); MAT_ASS_MAI(); MAT_ASS_BC() ; SOLVE(); OUTPUT_UCD() ; MAT_CO0: IU IALU 生成 MAT_CO: nde tem 生成 とりあえず から始まる節点番号を記憶
FEM3D 79 MAT_CO0: 全体構成 do cel= ICELTOT 節点相互の関係から ILUIALU を生成 (FID_ODE) enddo 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) /** ** MAT_CO0 **/ #nclude <stdo.h> #nclude "pfem_utl.h" #nclude "allocate.h etern FILE *fp_log; /*** eternal functons ***/ etern vod msort(nt* nt* nt); /*** statc functuons ***/ statc vod FID_TS_ODE (ntnt); vod MAT_CO0() { nt jkceln; nt nn2n3n4n5n6n7n; nt ; LU= 26; ILU=(KIT* )allocate_vector(seof(kit)); IALU=(KIT**)allocate_matr(seof(KIT)LU); for(=0;<;++) ILU[]=0; for(=0;<;++) for(j=0;j<lu;j++) IALU[][j]=0; LU: 各節点における非ゼロ非対角成分の最大数 ( 接続する節点数 ) 今の問題の場合はわかっているので このようにできる 不明の場合の実装 : レポート課題
FEM3D 行列コネクティビティ生成 : MAT_CO0(/4) /** ** MAT_CO0 **/ #nclude <stdo.h> #nclude "pfem_utl.h" #nclude "allocate.h etern FILE *fp_log; /*** eternal functons ***/ etern vod msort(nt* nt* nt); /*** statc functuons ***/ statc vod FID_TS_ODE (ntnt); vod MAT_CO0() { nt jkceln; nt nn2n3n4n5n6n7n; nt ; 変数名サイズ内容 ILU IALU [] []{[LU] 各節点の非零非対角成分数 各節点の非零非対角成分 ( 列番号 ) LU= 26; ILU=(KIT* )allocate_vector(seof(kit)); IALU=(KIT**)allocate_matr(seof(KIT)LU); for(=0;<;++) ILU[]=0; for(=0;<;++) for(j=0;j<lu;j++) IALU[][j]=0;
FEM3D 2 行列コネクティビティ生成 : MAT_CO0(2/4): から始まる番号 for( cel=0;cel< ICELTOT;cel++){ n=icelod[cel][0]; n2=icelod[cel][]; n3=icelod[cel][2]; n4=icelod[cel][3]; n5=icelod[cel][4]; n6=icelod[cel][5]; n7=icelod[cel][6]; n=icelod[cel][7]; FID_TS_ODE (nn2); FID_TS_ODE (nn3); FID_TS_ODE (nn4); FID_TS_ODE (nn5); FID_TS_ODE (nn6); FID_TS_ODE (nn7); FID_TS_ODE (nn); 7 5 6 4 3 2 FID_TS_ODE (n2n); FID_TS_ODE (n2n3); FID_TS_ODE (n2n4); FID_TS_ODE (n2n5); FID_TS_ODE (n2n6); FID_TS_ODE (n2n7); FID_TS_ODE (n2n); FID_TS_ODE (n3n); FID_TS_ODE (n3n2); FID_TS_ODE (n3n4); FID_TS_ODE (n3n5); FID_TS_ODE (n3n6); FID_TS_ODE (n3n7); FID_TS_ODE (n3n);
FEM3D 3 節点探索 :FID_TS_ODE ILUIAU 探索 : 一次元ではこの部分は手動 /*** *** FID_TS_ODE **/ statc vod FID_TS_ODE (nt pnt p2) { nt kkcou; for(kk=;kk<=ilu[p-];kk++){ f(p2 == IALU[p-][kk-]) return; cou=ilu[p-]+; IALU[p-][cou-]=p2; ILU[p-]=cou; return; 変数名サイズ内容 ILU IALU [] []{[LU] 各節点の非零非対角成分数 各節点の非零非対角成分 ( 列番号 )
FEM3D 4 節点探索 :FID_TS_ODE 一次元ではこの部分は手動 /*** *** FID_TS_ODE **/ statc vod FID_TS_ODE (nt pnt p2) { nt kkcou; for(kk=;kk<=ilu[p-];kk++){ f(p2 == IALU[p-][kk-]) return; 既に IALU に含まれている場合は 次のペアへ cou=ilu[p-]+; IALU[p-][cou-]=p2; ILU[p-]=cou; return; 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 一次元ではこの部分は手動 /*** *** FID_TS_ODE **/ statc vod FID_TS_ODE (nt pnt p2) { nt kkcou; for(kk=;kk<=ilu[p-];kk++){ f(p2 == IALU[p-][kk-]) return; cou=ilu[p-]+; IALU[p-][cou-]=p2; ILU[p-]=cou; return; 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) FID_TS_ODE (n4n); FID_TS_ODE (n4n2); FID_TS_ODE (n4n3); FID_TS_ODE (n4n5); FID_TS_ODE (n4n6); FID_TS_ODE (n4n7); FID_TS_ODE (n4n); FID_TS_ODE (n5n); FID_TS_ODE (n5n2); FID_TS_ODE (n5n3); FID_TS_ODE (n5n4); FID_TS_ODE (n5n6); FID_TS_ODE (n5n7); FID_TS_ODE (n5n); FID_TS_ODE (n6n); FID_TS_ODE (n6n2); FID_TS_ODE (n6n3); FID_TS_ODE (n6n4); FID_TS_ODE (n6n5); FID_TS_ODE (n6n7); FID_TS_ODE (n6n); 7 5 6 4 3 2 FID_TS_ODE (n7n); FID_TS_ODE (n7n2); FID_TS_ODE (n7n3); FID_TS_ODE (n7n4); FID_TS_ODE (n7n5); FID_TS_ODE (n7n6); FID_TS_ODE (n7n);
FEM3D 7 行列コネクティビティ生成 : MAT_CO0(4/4) FID_TS_ODE (nn); FID_TS_ODE (nn2); FID_TS_ODE (nn3); FID_TS_ODE (nn4); FID_TS_ODE (nn5); FID_TS_ODE (nn6); FID_TS_ODE (nn7); for(n=0;n<;n++){ =ILU[n]; for (k=0;k<;k++){ COL[k]=IALU[n][k]; msort(colcol2); for(k=;k>0;k--){ IALU[n][-k]= COL[COL2[k-]-]; 各節点において IALU[][k] が小さい番号から大きい番号に並ぶようにソート ( 単純なバブルソート ) せいぜい 00 程度のものをソートする
FEM3D CRS 形式への変換 :MAT_CO #nclude <stdo.h> #nclude "pfem_utl.h" #nclude "allocate.h" etern FILE* fp_log; vod MAT_CO() { nt kkk; ndelu=(kit*)allocate_vector(seof(kit)+); for(=0;<+;++) ndelu[]=0; for(=0;<;++){ ndelu[+]=ndelu[]+ilu[]; PLU=ndeLU[]; temlu=(kit*)allocate_vector(seof(kit)plu); for(=0;<;++){ for(k=0;k<ilu[];k++){ kk=k+ndelu[]; temlu[kk]=ialu[][k]-; deallocate_vector(ilu); deallocate_vector(ialu); C nde[ ] nde[0] 0 FORTRA nde(0) 0 k 0 ILU[ k] nde( ) ILU( k) k
FEM3D 9 CRS 形式への変換 :MAT_CO #nclude <stdo.h> #nclude "pfem_utl.h" #nclude "allocate.h" etern FILE* fp_log; vod MAT_CO() { nt kkk; ndelu=(kit*)allocate_vector(seof(kit)+); for(=0;<+;++) ndelu[]=0; for(=0;<;++){ ndelu[+]=ndelu[]+ilu[]; PLU=ndeLU[]; temlu=(kit*)allocate_vector(seof(kit)plu); for(=0;<;++){ for(k=0;k<ilu[];k++){ kk=k+ndelu[]; temlu[kk]=ialu[][k]-; PLU=nde[] tem のサイズ非ゼロ非対角成分総数 deallocate_vector(ilu); deallocate_vector(ialu); MAT_CO
FEM3D 90 CRS 形式への変換 :MAT_CO #nclude <stdo.h> #nclude "pfem_utl.h" #nclude "allocate.h" etern FILE* fp_log; vod MAT_CO() { nt kkk; ndelu=(kit*)allocate_vector(seof(kit)+); for(=0;<+;++) ndelu[]=0; for(=0;<;++){ ndelu[+]=ndelu[]+ilu[]; PLU=ndeLU[]; temlu=(kit*)allocate_vector(seof(kit)plu); for(=0;<;++){ for(k=0;k<ilu[];k++){ kk=k+ndelu[]; temlu[kk]=ialu[][k]-; tem に から始まる節点番号を記憶 deallocate_vector(ilu); deallocate_vector(ialu); MAT_CO
FEM3D 9 CRS 形式への変換 :MAT_CO #nclude <stdo.h> #nclude "pfem_utl.h" #nclude "allocate.h" etern FILE* fp_log; vod MAT_CO() { nt kkk; ndelu=(kit*)allocate_vector(seof(kit)+); for(=0;<+;++) ndelu[]=0; for(=0;<;++){ ndelu[+]=ndelu[]+ilu[]; PLU=ndeLU[]; temlu=(kit*)allocate_vector(seof(kit)plu); for(=0;<;++){ for(k=0;k<ilu[];k++){ kk=k+ndelu[]; temlu[kk]=ialu[][k]-; deallocate_vector(ilu); deallocate_vector(ialu); これらはもはや不要 MAT_CO
FEM3D 92 全体処理 /** program heat3d **/ #nclude <stdo.h> #nclude <stdlb.h> FILE* fp_log; #defne GLOBAL_VALUE_DEFIE #nclude "pfem_utl.h" //#nclude "solver.h" etern vod IPUT_CTL(); etern vod IPUT_GRID(); etern vod MAT_CO0(); etern vod MAT_CO(); etern vod MAT_ASS_MAI(); etern vod MAT_ASS_BC(); etern vod SOLVE(); etern vod OUTPUT_UCD(); nt man() { IPUT_CTL(); IPUT_GRID(); MAT_CO0(); MAT_CO(); MAT_ASS_MAI(); MAT_ASS_BC() ; SOLVE(); OUTPUT_UCD() ;
FEM3D 93 MAT_ASS_MAI: 全体構成 do kpn= 2 ガウス積分点番号 ( 方向 ) do jpn= 2 ガウス積分点番号 ( 方向 ) do pn= 2 ガウス積分点番号 ( 方向 ) ガウス積分点 ( 個 ) における形状関数 およびその 自然座標系 における微分の算出 enddo enddo enddo do cel= ICELTOT 要素ループ 節点の座標から ガウス積分点における 形状関数の 全体座標系 における微分 およびヤコビアンを算出 (JACOBI) do e= 局所節点番号 do je= 局所節点番号全体節点番号 :p jp A pjp のtemLUにおけるアドレス :kk j e do kpn= 2 ガウス積分点番号 ( 方向 ) do jpn= 2 ガウス積分点番号 ( 方向 ) do pn= 2 ガウス積分点番号 ( 方向 ) 要素積分 要素行列成分計算 全体行列への足しこみ enddo enddo enddo enddo enddo enddo e
FEM3D 94 係数行列 :MAT_ASS_MAI(/6) #nclude <stdo.h> #nclude <math.h> #nclude "pfem_utl.h" #nclude "allocate.h" etern FILE *fp_log; etern vod JACOBI(); vod MAT_ASS_MAI() { nt kkk; nt pjpkp; nt pnjpnkpn; nt cel; nt eje; nt SE; nt nn2n3n4n5n6n7n; double SH; double QPQMEPEMTPTM; double XX2X3X4X5X6X7X; double YY2Y3Y4Y5Y6Y7Y; double ZZ2Z3Z4Z5Z6Z7Z; double PXPYPZPXjPYjPZj; double COD0 QV0 QVC COEFj; double coef; KIT nodlocal[]; AMAT=(KREAL*) allocate_vector(seof(kreal)plu); 係数行列 ( 非零非対角成分 ) B =(KREAL*) allocate_vector(seof(kreal) ); 右辺ベクトル D =(KREAL*) allocate_vector(seof(kreal)); 解ベクトル X =(KREAL*) allocate_vector(seof(kreal)); 係数行列 ( 対角成分 ) for(=0;<plu;++) AMAT[]=0.0; for(=0;< ;++) B[]=0.0; for(=0;< ;++) D[]=0.0; for(=0;< ;++) X[]=0.0; WEI[0]=.0000000000e0; WEI[]=.0000000000e0; POS[0]= -0.5773502692e0; POS[]= 0.5773502692e0;
FEM3D 95 係数行列 :MAT_ASS_MAI(/6) #nclude <stdo.h> #nclude <math.h> #nclude "pfem_utl.h" #nclude "allocate.h" etern FILE *fp_log; etern vod JACOBI(); vod MAT_ASS_MAI() { nt kkk; nt pjpkp; nt pnjpnkpn; nt cel; nt eje; nt SE; nt nn2n3n4n5n6n7n; double SH; double QPQMEPEMTPTM; double XX2X3X4X5X6X7X; double YY2Y3Y4Y5Y6Y7Y; double ZZ2Z3Z4Z5Z6Z7Z; double PXPYPZPXjPYjPZj; double COD0 QV0 QVC COEFj; double coef; KIT nodlocal[]; AMAT=(KREAL*) allocate_vector(seof(kreal)plu); B =(KREAL*) allocate_vector(seof(kreal) ); D =(KREAL*) allocate_vector(seof(kreal)); X =(KREAL*) allocate_vector(seof(kreal)); for(=0;<plu;++) AMAT[]=0.0; for(=0;< ;++) B[]=0.0; for(=0;< ;++) D[]=0.0; for(=0;< ;++) X[]=0.0; WEI[0]=.0000000000e0; WEI[]=.0000000000e0; POS[0]= -0.5773502692e0; POS[]= 0.5773502692e0; POS: WEI: 積分点座標重み係数
FEM3D 96 係数行列 :MAT_ASS_MAI(2/6) /*** IIT. PQ - st-order dervatve of shape functon b QSI PE - st-order dervatve of shape functon b ETA PT - st-order dervatve of shape functon b ZET ***/ for(p=0;p<2;p++){ for(jp=0;jp<2;jp++){ for(kp=0;kp<2;kp++){ QP=.e0 + POS[p]; QM=.e0 - POS[p]; EP=.e0 + POS[jp]; EM=.e0 - POS[jp]; TP=.e0 + POS[kp]; TM=.e0 - POS[kp]; SHAPE[p][jp][kp][0]= Oth * QM * EM * TM; SHAPE[p][jp][kp][]= Oth * QP * EM * TM; SHAPE[p][jp][kp][2]= Oth * QP * EP * TM; SHAPE[p][jp][kp][3]= Oth * QM * EP * TM; SHAPE[p][jp][kp][4]= Oth * QM * EM * TP; SHAPE[p][jp][kp][5]= Oth * QP * EM * TP; SHAPE[p][jp][kp][6]= Oth * QP * EP * TP; SHAPE[p][jp][kp][7]= Oth * QM * EP * TP;
FEM3D 97 係数行列 :MAT_ASS_MAI(2/6) /*** IIT. PQ - st-order dervatve of shape functon b QSI PE - st-order dervatve of shape functon b ETA PT - st-order dervatve of shape functon b ZET ***/ for(p=0;p<2;p++){ for(jp=0;jp<2;jp++){ for(kp=0;kp<2;kp++){ QP=.e0 + POS[p]; QM=.e0 - POS[p]; EP=.e0 + POS[jp]; EM=.e0 - POS[jp]; TP=.e0 + POS[kp]; TM=.e0 - POS[kp]; SHAPE[p][jp][kp][0]= Oth * QM * EM * TM; SHAPE[p][jp][kp][]= Oth * QP * EM * TM; SHAPE[p][jp][kp][2]= Oth * QP * EP * TM; SHAPE[p][jp][kp][3]= Oth * QM * EP * TM; SHAPE[p][jp][kp][4]= Oth * QM * EM * TP; SHAPE[p][jp][kp][5]= Oth * QP * EM * TP; SHAPE[p][jp][kp][6]= Oth * QP * EP * TP; SHAPE[p][jp][kp][7]= Oth * QM * EP * TP; QP EP TP k QM j EM j j TMk k k
FEM3D 9 係数行列 :MAT_ASS_MAI(2/6) /*** IIT. PQ - st-order dervatve of shape functon b QSI PE - st-order dervatve of shape functon b ETA PT - st-order dervatve of shape functon b ZET ***/ for(p=0;p<2;p++){ for(jp=0;jp<2;jp++){ for(kp=0;kp<2;kp++){ QP=.e0 + POS[p]; QM=.e0 - POS[p]; EP=.e0 + POS[jp]; EM=.e0 - POS[jp]; TP=.e0 + POS[kp]; TM=.e0 - POS[kp]; SHAPE[p][jp][kp][0]= Oth * QM * EM * TM; SHAPE[p][jp][kp][]= Oth * QP * EM * TM; SHAPE[p][jp][kp][2]= Oth * QP * EP * TM; SHAPE[p][jp][kp][3]= Oth * QM * EP * TM; SHAPE[p][jp][kp][4]= Oth * QM * EM * TP; SHAPE[p][jp][kp][5]= Oth * QP * EM * TP; SHAPE[p][jp][kp][6]= Oth * QP * EP * TP; SHAPE[p][jp][kp][7]= Oth * QM * EP * TP; 7 5 6 4 2 3
FEM3D 99 係数行列 :MAT_ASS_MAI(2/6) /*** IIT. PQ - st-order dervatve of shape functon b QSI PE - st-order dervatve of shape functon b ETA PT - st-order dervatve of shape functon b ZET ***/ for(p=0;p<2;p++){ for(jp=0;jp<2;jp++){ for(kp=0;kp<2;kp++){ QP=.e0 + POS[p]; QM=.e0 - POS[p]; EP=.e0 + POS[jp]; EM=.e0 - POS[jp]; TP=.e0 + POS[kp]; TM=.e0 - POS[kp]; SHAPE[p][jp][kp][0]= Oth * QM * EM * TM; SHAPE[p][jp][kp][]= Oth * QP * EM * TM; SHAPE[p][jp][kp][2]= Oth * QP * EP * TM; SHAPE[p][jp][kp][3]= Oth * QM * EP * TM; SHAPE[p][jp][kp][4]= Oth * QM * EM * TP; SHAPE[p][jp][kp][5]= Oth * QP * EM * TP; SHAPE[p][jp][kp][6]= Oth * QP * EP * TP; SHAPE[p][jp][kp][7]= Oth * QM * EP * TP; ) ( ) ( ) ( ) ( 4 3 2 ) ( ) ( ) ( ) ( 7 6 5
FEM3D 00 係数行列 :MAT_ASS_MAI(3/6) PQ[jp][kp][0]= - Oth * EM * TM; PQ[jp][kp][]= + Oth * EM * TM; PQ[jp][kp][2]= + Oth * EP * TM; PQ[jp][kp][3]= - Oth * EP * TM; PQ[jp][kp][4]= - Oth * EM * TP; PQ[jp][kp][5]= + Oth * EM * TP; PQ[jp][kp][6]= + Oth * EP * TP; PQ[jp][kp][7]= - Oth * EP * TP; PE[p][kp][0]= - Oth * QM * TM; PE[p][kp][]= - Oth * QP * TM; PE[p][kp][2]= + Oth * QP * TM; PE[p][kp][3]= + Oth * QM * TM; PE[p][kp][4]= - Oth * QM * TP; PE[p][kp][5]= - Oth * QP * TP; PE[p][kp][6]= + Oth * QP * TP; PE[p][kp][7]= + Oth * QM * TP; PT[p][jp][0]= - Oth * QM * EM; PT[p][jp][]= - Oth * QP * EM; PT[p][jp][2]= - Oth * QP * EP; PT[p][jp][3]= - Oth * QM * EP; PT[p][jp][4]= + Oth * QM * EM; PT[p][jp][5]= + Oth * QP * EM; PT[p][jp][6]= + Oth * QP * EP; PT[p][jp][7]= + Oth * QM * EP; for( cel=0;cel< ICELTOT;cel++){ COD0= COD; n=icelod[cel][0]; n2=icelod[cel][]; n3=icelod[cel][2]; n4=icelod[cel][3]; n5=icelod[cel][4]; n6=icelod[cel][5]; n7=icelod[cel][6]; n=icelod[cel][7]; ( 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[jp][kp][0]= - Oth * EM * TM; PQ[jp][kp][]= + Oth * EM * TM; PQ[jp][kp][2]= + Oth * EP * TM; PQ[jp][kp][3]= - Oth * EP * TM; PQ[jp][kp][4]= - Oth * EM * TP; PQ[jp][kp][5]= + Oth * EM * TP; PQ[jp][kp][6]= + Oth * EP * TP; PQ[jp][kp][7]= - Oth * EP * TP; PE[p][kp][0]= - Oth * QM * TM; PE[p][kp][]= - Oth * QP * TM; PE[p][kp][2]= + Oth * QP * TM; PE[p][kp][3]= + Oth * QM * TM; PE[p][kp][4]= - Oth * QM * TP; PE[p][kp][5]= - Oth * QP * TP; PE[p][kp][6]= + Oth * QP * TP; PE[p][kp][7]= + Oth * QM * TP; PT[p][jp][0]= - Oth * QM * EM; PT[p][jp][]= - Oth * QP * EM; PT[p][jp][2]= - Oth * QP * EP; PT[p][jp][3]= - Oth * QM * EP; PT[p][jp][4]= + Oth * QM * EM; PT[p][jp][5]= + Oth * QP * EM; PT[p][jp][6]= + Oth * QP * EP; PT[p][jp][7]= + Oth * QM * EP; for( cel=0;cel< ICELTOT;cel++){ COD0= COD; n=icelod[cel][0]; n2=icelod[cel][]; n3=icelod[cel][2]; n4=icelod[cel][3]; n5=icelod[cel][4]; n6=icelod[cel][5]; n7=icelod[cel][6]; n=icelod[cel][7]; 7 5 6 4 3 2
FEM3D 02 係数行列 :MAT_ASS_MAI(4/6) nodlocal[0]= n; nodlocal[]= n2; nodlocal[2]= n3; nodlocal[3]= n4; nodlocal[4]= n5; nodlocal[5]= n6; nodlocal[6]= n7; nodlocal[7]= n; X=XYZ[n-][0]; X2=XYZ[n2-][0]; X3=XYZ[n3-][0]; X4=XYZ[n4-][0]; X5=XYZ[n5-][0]; X6=XYZ[n6-][0]; X7=XYZ[n7-][0]; X=XYZ[n-][0]; Y=XYZ[n-][]; Y2=XYZ[n2-][]; Y3=XYZ[n3-][]; Y4=XYZ[n4-][]; Y5=XYZ[n5-][]; Y6=XYZ[n6-][]; Y7=XYZ[n7-][]; Y=XYZ[n-][]; 節点の節点番号 7 5 6 4 2 3 QVC= Oth*(X+X2+X3+X4+X5+X6+X7+X+Y+Y2+Y3+Y4+Y5+Y6+Y7+Y); Z=XYZ[n-][2]; Z2=XYZ[n2-][2]; Z3=XYZ[n3-][2]; Z4=XYZ[n4-][2]; Z5=XYZ[n5-][2]; Z6=XYZ[n6-][2]; Z7=XYZ[n7-][2]; Z=XYZ[n-][2]; 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[0]= n; nodlocal[]= n2; nodlocal[2]= n3; nodlocal[3]= n4; nodlocal[4]= n5; nodlocal[5]= n6; nodlocal[6]= n7; nodlocal[7]= n; X=XYZ[n-][0]; X2=XYZ[n2-][0]; X3=XYZ[n3-][0]; X4=XYZ[n4-][0]; X5=XYZ[n5-][0]; X6=XYZ[n6-][0]; X7=XYZ[n7-][0]; X=XYZ[n-][0]; Y=XYZ[n-][]; Y2=XYZ[n2-][]; Y3=XYZ[n3-][]; Y4=XYZ[n4-][]; Y5=XYZ[n5-][]; Y6=XYZ[n6-][]; Y7=XYZ[n7-][]; Y=XYZ[n-][]; 節点の X 座標 節点の Y 座標 2 5 6 4 7 3 座標値 : 節点番号から 引く QVC= Oth*(X+X2+X3+X4+X5+X6+X7+X+Y+Y2+Y3+Y4+Y5+Y6+Y7+Y); Z=XYZ[n-][2]; Z2=XYZ[n2-][2]; Z3=XYZ[n3-][2]; Z4=XYZ[n4-][2]; Z5=XYZ[n5-][2]; Z6=XYZ[n6-][2]; Z7=XYZ[n7-][2]; Z=XYZ[n-][2]; 節点の Z 座標 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[0]= n; nodlocal[]= n2; nodlocal[2]= n3; nodlocal[3]= n4; nodlocal[4]= n5; nodlocal[5]= n6; nodlocal[6]= n7; nodlocal[7]= n; X=XYZ[n-][0]; X2=XYZ[n2-][0]; X3=XYZ[n3-][0]; X4=XYZ[n4-][0]; X5=XYZ[n5-][0]; X6=XYZ[n6-][0]; X7=XYZ[n7-][0]; X=XYZ[n-][0]; Y=XYZ[n-][]; Y2=XYZ[n2-][]; Y3=XYZ[n3-][]; Y4=XYZ[n4-][]; Y5=XYZ[n5-][]; Y6=XYZ[n6-][]; Y7=XYZ[n7-][]; Y=XYZ[n-][]; 節点の X 座標 節点の Y 座標 2 5 6 4 7 3 座標値 : 節点番号から 引く QVC= Oth*(X+X2+X3+X4+X5+X6+X7+X+Y+Y2+Y3+Y4+Y5+Y6+Y7+Y); Z=XYZ[n-][2]; Z2=XYZ[n2-][2]; Z3=XYZ[n3-][2]; Z4=XYZ[n4-][2]; Z5=XYZ[n5-][2]; Z6=XYZ[n6-][2]; Z7=XYZ[n7-][2]; Z=XYZ[n-][2]; T T 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); Q QVOL C C T Q 0 体積当たり発熱量は位置 ( メッシュの中心の座標 c c ) に依存
FEM3D 05 係数行列 :MAT_ASS_MAI(4/6) nodlocal[0]= n; nodlocal[]= n2; nodlocal[2]= n3; nodlocal[3]= n4; nodlocal[4]= n5; nodlocal[5]= n6; nodlocal[6]= n7; nodlocal[7]= n; X=XYZ[n-][0]; X2=XYZ[n2-][0]; X3=XYZ[n3-][0]; X4=XYZ[n4-][0]; X5=XYZ[n5-][0]; X6=XYZ[n6-][0]; X7=XYZ[n7-][0]; X=XYZ[n-][0]; Y=XYZ[n-][]; Y2=XYZ[n2-][]; Y3=XYZ[n3-][]; Y4=XYZ[n4-][]; Y5=XYZ[n5-][]; Y6=XYZ[n6-][]; Y7=XYZ[n7-][]; Y=XYZ[n-][]; 2 5 6 4 7 3 座標値 : 節点番号から 引く QVC= Oth*(X+X2+X3+X4+X5+X6+X7+X+Y+Y2+Y3+Y4+Y5+Y6+Y7+Y); Z=XYZ[n-][2]; Z2=XYZ[n2-][2]; Z3=XYZ[n3-][2]; Z4=XYZ[n4-][2]; Z5=XYZ[n5-][2]; Z6=XYZ[n6-][2]; Z7=XYZ[n7-][2]; Z=XYZ[n-][2]; T QVC C C 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); Q T QVOL C C T Q 0
FEM3D 06 係数行列 :MAT_ASS_MAI(4/6) nodlocal[0]= n; nodlocal[]= n2; nodlocal[2]= n3; nodlocal[3]= n4; nodlocal[4]= n5; nodlocal[5]= n6; nodlocal[6]= n7; nodlocal[7]= n; X=XYZ[n-][0]; X2=XYZ[n2-][0]; X3=XYZ[n3-][0]; X4=XYZ[n4-][0]; X5=XYZ[n5-][0]; X6=XYZ[n6-][0]; X7=XYZ[n7-][0]; X=XYZ[n-][0]; Y=XYZ[n-][]; Y2=XYZ[n2-][]; Y3=XYZ[n3-][]; Y4=XYZ[n4-][]; Y5=XYZ[n5-][]; Y6=XYZ[n6-][]; Y7=XYZ[n7-][]; Y=XYZ[n-][]; QVC= Oth*(X+X2+X3+X4+X5+X6+X7+X+Y+Y2+Y3+Y4+Y5+Y6+Y7+Y); Z=XYZ[n-][2]; Z2=XYZ[n2-][2]; Z3=XYZ[n3-][2]; Z4=XYZ[n4-][2]; Z5=XYZ[n5-][2]; Z6=XYZ[n6-][2]; Z7=XYZ[n7-][2]; Z=XYZ[n-][2]; 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) #nclude <stdo.h> #nclude <math.h> #nclude "precson.h" #nclude "allocate.h" /*** *** JACOBI ***/ vod JACOBI( KREAL DETJ[2][2][2] KREAL PQ[2][2][]KREAL PE[2][2][]KREAL PT[2][2][] KREAL PX[2][2][2][]KREAL PY[2][2][2][]KREAL PZ[2][2][2][] KREAL XKREAL X2KREAL X3KREAL X4KREAL X5KREAL X6KREAL X7KREAL X KREAL YKREAL Y2KREAL Y3KREAL Y4KREAL Y5KREAL Y6KREAL Y7KREAL Y KREAL ZKREAL Z2KREAL Z3KREAL Z4KREAL Z5KREAL Z6KREAL Z7KREAL Z) { /** calculates JACOBIA & IVERSE JACOBIA d/d d/d & d/d **/ nt pjpkp; double dxdqdydqdzdqdxdedydedzdedxdtdydtdzdt; double coef; double aa2a3a2a22a23a3a32a33; for(p=0;p<2;p++){ for(jp=0;jp<2;jp++){ for(kp=0;kp<2;kp++){ PX[p][jp][kp][0]=0.0; PX[p][jp][kp][]=0.0; PX[p][jp][kp][2]=0.0; PX[p][jp][kp][3]=0.0; PX[p][jp][kp][4]=0.0; PX[p][jp][kp][5]=0.0; PX[p][jp][kp][6]=0.0; PX[p][jp][kp][7]=0.0; l l l 入力 l ~ 出力 l l l det J l l l 各ガウス積分点 [p][jp][kp] における値
FEM3D 0 自然座標系における偏微分 (/4) 偏微分の公式より以下のようになる : ) ( ) ( ) ( は定義より簡単に求められるがを実際の計算で使用する
FEM3D 09 自然座標系における偏微分 (2/4) マトリックス表示すると : J 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) の定義より簡単に求められる J J J J J J J J J 33 32 3 23 22 2 3 2
FEM3D JACOBI(2/4) PY[p][jp][kp][0]=0.0; PY[p][jp][kp][]=0.0; PY[p][jp][kp][2]=0.0; PY[p][jp][kp][3]=0.0; PY[p][jp][kp][4]=0.0; PY[p][jp][kp][5]=0.0; PY[p][jp][kp][6]=0.0; PY[p][jp][kp][7]=0.0; PZ[p][jp][kp][0]=0.0; PZ[p][jp][kp][]=0.0; PZ[p][jp][kp][2]=0.0; PZ[p][jp][kp][3]=0.0; PZ[p][jp][kp][4]=0.0; PZ[p][jp][kp][5]=0.0; PZ[p][jp][kp][6]=0.0; PZ[p][jp][kp][7]=0.0; J J J J 2 3 J J J 2 22 32 J J J 3 23 33 /** **/ DETERMIAT of the JACOBIA dxdq = PQ[jp][kp][0]*X + PQ[jp][kp][]*X2 + PQ[jp][kp][2]*X3 + PQ[jp][kp][3]*X4 + PQ[jp][kp][4]*X5 + PQ[jp][kp][5]*X6 + PQ[jp][kp][6]*X7 + PQ[jp][kp][7]*X; dydq = PQ[jp][kp][0]*Y + PQ[jp][kp][]*Y2 + PQ[jp][kp][2]*Y3 + PQ[jp][kp][3]*Y4 + PQ[jp][kp][4]*Y5 + PQ[jp][kp][5]*Y6 + PQ[jp][kp][6]*Y7 + PQ[jp][kp][7]*Y; dzdq = PQ[jp][kp][0]*Z + PQ[jp][kp][]*Z2 + PQ[jp][kp][2]*Z3 + PQ[jp][kp][3]*Z4 + PQ[jp][kp][4]*Z5 + PQ[jp][kp][5]*Z6 + PQ[jp][kp][6]*Z7 + PQ[jp][kp][7]*Z; dxde = PE[p][kp][0]*X + PE[p][kp][]*X2 + PE[p][kp][2]*X3 + PE[p][kp][3]*X4 + PE[p][kp][4]*X5 + PE[p][kp][5]*X6 + PE[p][kp][6]*X7 + PE[p][kp][7]*X; dxdq J dydq J dzdq J 2 3
FEM3D 2 /** IVERSE JACOBIA **/ JACOBI(3/4) dyde = PE[p][kp][0]*Y + PE[p][kp][]*Y2 + PE[p][kp][2]*Y3 + PE[p][kp][3]*Y4 + PE[p][kp][4]*Y5 + PE[p][kp][5]*Y6 + PE[p][kp][6]*Y7 + PE[p][kp][7]*Y; dzde = PE[p][kp][0]*Z + PE[p][kp][]*Z2 + PE[p][kp][2]*Z3 + PE[p][kp][3]*Z4 + PE[p][kp][4]*Z5 + PE[p][kp][5]*Z6 + PE[p][kp][6]*Z7 + PE[p][kp][7]*Z; dxdt = PT[p][jp][0]*X + PT[p][jp][]*X2 + PT[p][jp][2]*X3 + PT[p][jp][3]*X4 + PT[p][jp][4]*X5 + PT[p][jp][5]*X6 + PT[p][jp][6]*X7 + PT[p][jp][7]*X; dydt = PT[p][jp][0]*Y + PT[p][jp][]*Y2 + PT[p][jp][2]*Y3 + PT[p][jp][3]*Y4 + PT[p][jp][4]*Y5 + PT[p][jp][5]*Y6 + PT[p][jp][6]*Y7 + PT[p][jp][7]*Y; dzdt = PT[p][jp][0]*Z + PT[p][jp][]*Z2 + PT[p][jp][2]*Z3 + PT[p][jp][3]*Z4 + PT[p][jp][4]*Z5 + PT[p][jp][5]*Z6 + PT[p][jp][6]*Z7 + PT[p][jp][7]*Z; DETJ[p][jp][kp]= dxdq*(dyde*dzdt-dzde*dydt) + dydq*(dzde*dxdt-dxde*dzdt) + dzdq*(dxde*dydt-dyde*dxdt); coef=.0 / DETJ[p][jp][kp]; 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 ); DETJ[p][jp][kp]=fabs(DETJ[p][jp][kp]); J J J J 2 3 J J J 2 22 32 J J J 3 23 33
FEM3D 3 自然座標系における偏微分 (4/4) 従って下記のように偏微分を計算できる ヤコビアン (3 3 行列 ) の逆行列を求める J
FEM3D 4 /** IVERSE JACOBIA JACOBI(3/4) dyde = PE[p][kp][0]*Y + PE[p][kp][]*Y2 + PE[p][kp][2]*Y3 + PE[p][kp][3]*Y4 + PE[p][kp][4]*Y5 + PE[p][kp][5]*Y6 + PE[p][kp][6]*Y7 + PE[p][kp][7]*Y; dzde = PE[p][kp][0]*Z + PE[p][kp][]*Z2 + PE[p][kp][2]*Z3 + PE[p][kp][3]*Z4 + PE[p][kp][4]*Z5 + PE[p][kp][5]*Z6 + PE[p][kp][6]*Z7 + PE[p][kp][7]*Z; dxdt = PT[p][jp][0]*X + PT[p][jp][]*X2 + PT[p][jp][2]*X3 + PT[p][jp][3]*X4 + PT[p][jp][4]*X5 + PT[p][jp][5]*X6 + PT[p][jp][6]*X7 + PT[p][jp][7]*X; dydt = PT[p][jp][0]*Y + PT[p][jp][]*Y2 + PT[p][jp][2]*Y3 + PT[p][jp][3]*Y4 + PT[p][jp][4]*Y5 + PT[p][jp][5]*Y6 + PT[p][jp][6]*Y7 + PT[p][jp][7]*Y; dzdt = PT[p][jp][0]*Z + PT[p][jp][]*Z2 + PT[p][jp][2]*Z3 + PT[p][jp][3]*Z4 + PT[p][jp][4]*Z5 + PT[p][jp][5]*Z6 + PT[p][jp][6]*Z7 + PT[p][jp][7]*Z; DETJ[p][jp][kp]= dxdq*(dyde*dzdt-dzde*dydt) + dydq*(dzde*dxdt-dxde*dzdt) + dzdq*(dxde*dydt-dyde*dxdt); **/ coef=.0 / DETJ[p][jp][kp]; 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[p][jp][kp]=fabs(DETJ[p][jp][kp]);
FEM3D 5 JACOBI(4/4) /** set the d/dx d/dy & d/dz components **/ PX[p][jp][kp][0]=a*PQ[jp][kp][0]+a2*PE[p][kp][0]+a3*PT[p][jp][0]; PX[p][jp][kp][]=a*PQ[jp][kp][]+a2*PE[p][kp][]+a3*PT[p][jp][]; PX[p][jp][kp][2]=a*PQ[jp][kp][2]+a2*PE[p][kp][2]+a3*PT[p][jp][2]; PX[p][jp][kp][3]=a*PQ[jp][kp][3]+a2*PE[p][kp][3]+a3*PT[p][jp][3]; PX[p][jp][kp][4]=a*PQ[jp][kp][4]+a2*PE[p][kp][4]+a3*PT[p][jp][4]; PX[p][jp][kp][5]=a*PQ[jp][kp][5]+a2*PE[p][kp][5]+a3*PT[p][jp][5]; PX[p][jp][kp][6]=a*PQ[jp][kp][6]+a2*PE[p][kp][6]+a3*PT[p][jp][6]; PX[p][jp][kp][7]=a*PQ[jp][kp][7]+a2*PE[p][kp][7]+a3*PT[p][jp][7]; PY[p][jp][kp][0]=a2*PQ[jp][kp][0]+a22*PE[p][kp][0]+a23*PT[p][jp][0]; PY[p][jp][kp][]=a2*PQ[jp][kp][]+a22*PE[p][kp][]+a23*PT[p][jp][]; PY[p][jp][kp][2]=a2*PQ[jp][kp][2]+a22*PE[p][kp][2]+a23*PT[p][jp][2]; PY[p][jp][kp][3]=a2*PQ[jp][kp][3]+a22*PE[p][kp][3]+a23*PT[p][jp][3]; PY[p][jp][kp][4]=a2*PQ[jp][kp][4]+a22*PE[p][kp][4]+a23*PT[p][jp][4]; PY[p][jp][kp][5]=a2*PQ[jp][kp][5]+a22*PE[p][kp][5]+a23*PT[p][jp][5]; PY[p][jp][kp][6]=a2*PQ[jp][kp][6]+a22*PE[p][kp][6]+a23*PT[p][jp][6]; PY[p][jp][kp][7]=a2*PQ[jp][kp][7]+a22*PE[p][kp][7]+a23*PT[p][jp][7]; PZ[p][jp][kp][0]=a3*PQ[jp][kp][0]+a32*PE[p][kp][0]+a33*PT[p][jp][0]; PZ[p][jp][kp][]=a3*PQ[jp][kp][]+a32*PE[p][kp][]+a33*PT[p][jp][]; PZ[p][jp][kp][2]=a3*PQ[jp][kp][2]+a32*PE[p][kp][2]+a33*PT[p][jp][2]; PZ[p][jp][kp][3]=a3*PQ[jp][kp][3]+a32*PE[p][kp][3]+a33*PT[p][jp][3]; PZ[p][jp][kp][4]=a3*PQ[jp][kp][4]+a32*PE[p][kp][4]+a33*PT[p][jp][4]; PZ[p][jp][kp][5]=a3*PQ[jp][kp][5]+a32*PE[p][kp][5]+a33*PT[p][jp][5]; PZ[p][jp][kp][6]=a3*PQ[jp][kp][6]+a32*PE[p][kp][6]+a33*PT[p][jp][6]; PZ[p][jp][kp][7]=a3*PQ[jp][kp][7]+a32*PE[p][kp][7]+a33*PT[p][jp][7]; a a a a a a a a a 33 32 3 23 22 2 3 2
FEM3D 6 係数行列 :MAT_ASS_MAI(5/6) /** COSTRUCT the GLOBAL MATRIX **/ for(e=0;e<;e++){ p=nodlocal[e]; for(je=0;je<;je++){ jp=nodlocal[je]; kk=-; f( jp!= p ){ S=ndeLU[p-]; E=ndeLU[p ]; for( k=s;k<e;k++){ f( temlu[k] == jp- ){ kk=k; break; 全体行列の非対角成分 A p jp kk:tem におけるアドレス p= nodlocal[e] jp= nodlocal[je] 7 5 6 4 3 から始まる節点番号 2
FEM3D 7 要素マトリクス : 行列 j j k j 7 5 6 4 3 2
FEM3D 係数行列 :MAT_ASS_MAI(5/6) /** COSTRUCT the GLOBAL MATRIX **/ for(e=0;e<;e++){ p=nodlocal[e]; for(je=0;je<;je++){ jp=nodlocal[je]; kk=-; f( jp!= p ){ S=ndeLU[p-]; E=ndeLU[p ]; for( k=s;k<e;k++){ f( temlu[k] == jp- ){ kk=k; break; 要素マトリクス ( e ~j e ) 全体マトリクス ( p ~j p ) の関係 kk:temlu におけるアドレス 7 5 6 4 3 2
FEM3D 9 係数行列 :MAT_ASS_MAI(6/6) QV0= 0.e0; COEFj= 0.e0; for(kpn=0;kpn<2;kpn++){ for(jpn=0;jpn<2;jpn++){ for(pn=0;pn<2;pn++){ coef= fabs(detj[pn][jpn][kpn])*wei[pn]*wei[jpn]*wei[kpn]; PX= PX[pn][jpn][kpn][e]; PY= PY[pn][jpn][kpn][e]; PZ= PZ[pn][jpn][kpn][e]; PXj= PX[pn][jpn][kpn][je]; PYj= PY[pn][jpn][kpn][je]; PZj= PZ[pn][jpn][kpn][je]; COEFj+= coef*cod0*(px*pxj+py*pyj+pz*pzj); SH= SHAPE[pn][jpn][kpn][e]; QV0+= SH * QVOL * coef; f (jp==p) { D[p-]+= COEFj; B[p-]+= QV0*QVC; f (jp!= p) { AMAT[kk]+= COEFj; j j j det J ddd
FEM3D 20 係数行列 :MAT_ASS_MAI(6/6) QV0= 0.e0; COEFj= 0.e0; for(kpn=0;kpn<2;kpn++){ for(jpn=0;jpn<2;jpn++){ for(pn=0;pn<2;pn++){ coef= fabs(detj[pn][jpn][kpn])*wei[pn]*wei[jpn]*wei[kpn]; PX= PX[pn][jpn][kpn][e]; PY= PY[pn][jpn][kpn][e]; PZ= PZ[pn][jpn][kpn][e]; PXj= PX[pn][jpn][kpn][je]; PYj= PY[pn][jpn][kpn][je]; PZj= PZ[pn][jpn][kpn][je]; COEFj+= coef*cod0*(px*pxj+py*pyj+pz*pzj); SH= SHAPE[pn][jpn][kpn][e]; QV0+= SH * QVOL * coef; f (jp==p) { D[p-]+= COEFj; B[p-]+= QV0*QVC; f (jp!= p) { AMAT[kk]+= COEFj; d d d J j j j det k k j k j M j L f W W W d d d f I ) ( ) (
FEM3D 2 係数行列 :MAT_ASS_MAI(6/6) QV0= 0.e0; COEFj= 0.e0; for(kpn=0;kpn<2;kpn++){ for(jpn=0;jpn<2;jpn++){ for(pn=0;pn<2;pn++){ coef= fabs(detj[pn][jpn][kpn])*wei[pn]*wei[jpn]*wei[kpn]; PX= PX[pn][jpn][kpn][e]; PY= PY[pn][jpn][kpn][e]; PZ= PZ[pn][jpn][kpn][e]; PXj= PX[pn][jpn][kpn][je]; PYj= PY[pn][jpn][kpn][je]; PZj= PZ[pn][jpn][kpn][je]; COEFj+= coef*cod0*(px*pxj+py*pyj+pz*pzj); SH= SHAPE[pn][jpn][kpn][e]; QV0+= SH * QVOL * coef; f (jp==p) { D[p-]+= COEFj; B[p-]+= QV0*QVC; f (jp!= p) { AMAT[kk]+= COEFj; d d d J j j j det k k j k j M j L f W W W d d d f I ) ( ) ( k j k j J W W W coef det
FEM3D 22 係数行列 :MAT_ASS_MAI(6/6) QV0= 0.e0; COEFj= 0.e0; for(kpn=0;kpn<2;kpn++){ for(jpn=0;jpn<2;jpn++){ for(pn=0;pn<2;pn++){ coef= fabs(detj[pn][jpn][kpn])*wei[pn]*wei[jpn]*wei[kpn]; PX= PX[pn][jpn][kpn][e]; PY= PY[pn][jpn][kpn][e]; PZ= PZ[pn][jpn][kpn][e]; PXj= PX[pn][jpn][kpn][je]; PYj= PY[pn][jpn][kpn][je]; PZj= PZ[pn][jpn][kpn][je]; COEFj+= coef*cod0*(px*pxj+py*pyj+pz*pzj); SH= SHAPE[pn][jpn][kpn][e]; QV0+= SH * QVOL * coef; j k j j f (jp==p) { D[p-]+= COEFj; B[p-]+= QV0*QVC; f (jp!= p) { AMAT[kk]+= COEFj;
FEM3D 23 係数行列 :MAT_ASS_MAI(6/6) QV0= 0.e0; COEFj= 0.e0; for(kpn=0;kpn<2;kpn++){ for(jpn=0;jpn<2;jpn++){ for(pn=0;pn<2;pn++){ coef= fabs(detj[pn][jpn][kpn])*wei[pn]*wei[jpn]*wei[kpn]; PX= PX[pn][jpn][kpn][e]; PY= PY[pn][jpn][kpn][e]; k ( e) ( e) f ( e PZ= PZ[pn][jpn][kpn][e]; PXj= PX[pn][jpn][kpn][je]; PYj= PY[pn][jpn][kpn][je]; PZj= PZ[pn][jpn][kpn][je]; COEFj+= coef*cod0*(px*pxj+py*pyj+pz*pzj); SH= SHAPE[pn][jpn][kpn][e]; QV0+= SH * QVOL * coef; f (jp==p) { D[p-]+= COEFj; B[p-]+= QV0*QVC; f (jp!= p) { AMAT[kk]+= COEFj; ) ( e) T f Q dv Q V QVOL C C QVC C C QV 0 QVOL f ( e) V QV 0QVC T dv
FEM3D 24 MAT_ASS_BC: 全体構成 do = 節点ループ ( ディリクレ ) 境界条件を設定する節点をマーク (IWKX) enddo do = 節点ループ f (IWKX().eq.) then マークされた節点だったら対応する右辺ベクトル (B) の成分 対角成分 (D) の成分の修正 ( 行 列 ) do k= nde(-)+ nde() 対応する非零非対角成分 (AMAT) の成分の修正 ( 行 ) enddo endf enddo Z T=0@Z= ma do = 節点ループ do k= ndelu (-)+ nde () f (IWKX(tem (k)).eq.) then 対応する非零非対角成分の節点がマークされていたら対応する右辺ベクトル 非零非対角成分 (AMAT) の成分の修正 ( 列 ) endf enddo enddo X Y X Z Y
FEM3D 25 境界条件 :MAT_ASS_BC(/2) #nclude <stdo.h> #nclude <stdlb.h> #nclude <strng.h> #nclude "pfem_utl.h" #nclude "allocate.h" etern FILE *fp_log; vod MAT_ASS_BC() { nt jknbb0cel; nt nn2n3n4n5n6n7n; nt qq2q3q4q5q6q7q; nt SE; double STRESSVAL; IWKX=(KIT**) allocate_matr(seof(kit)2); for(=0;<;++) for(j=0;j<2;j++) IWKX[][j]=0; /** Z=Zma **/ for(n=0;n<;n++) IWKX[n][0]=0; b0=-; 節点グループ名が Zma である節点 n( から始まる ) において : IWKX[n-][0]= とする for( b0=0;b0<odgrptot;b0++){ f( strcmp(odgrp_ame[b0].name"zma") == 0 ) break; for( b=odgrp_idex[b0];b<odgrp_idex[b0+];b++){ n=odgrp_item[b]; IWKX[n-][0]=;
FEM3D 26 境界条件 :MAT_ASS_BC(2/2) for(n=0;n<;n++){ f( IWKX[n][0] == ){ B[n]= 0.e0; D[n]=.e0; for(k=ndelu[n];k<ndelu[n+];k++){ AMAT[k]= 0.e0; for(n=0;n<;n++){ for(k=ndelu[n];k<ndelu[n+];k++){ f (IWKX[temLU[k]][0] == ) { AMAT[k]= 0.e0;
FEM3D 27 対象とする問題 : 一次元熱伝導問題 体積当たり一様発熱 Q T Q 0 =0 ( mn ) = ma 一様な : 断面積 A 熱伝導率 体積当たり一様発熱 ( 時間当たり ) QL -3 T - 境界条件 =0 :=0 ( 固定 ) T = ma : 0( 断熱 ) Q 復習 : 一次元問題
FEMD 2 =0 で成立する方程式 T 0 =0 体積当たり一様発熱 Q T Q 0 =0 ( mn ) = ma 一様な : 断面積 A 熱伝導率 体積当たり一様発熱 ( 時間当たり ) QL -3 T - 境界条件 =0 :=0 ( 固定 ) T = ma : 0( 断熱 ) Q 復習 : 一次元問題
FEMD 29 プログラム :d.c(6/6) 第一種境界条件 @=0 /* // +---------------------+ // BOUDARY condtons // +---------------------+ */ /* X=Xmn */ =0; js= Inde[]; AMat[jS]= 0.0; Dag[ ]=.0; Rhs [ ]= 0.0; T 0 =0 対角成分 = 右辺 =0 非対角成分 =0 for(k=0;k<plu;k++){ f(item[k]==0){amat[k]=0.0; 復習 : 一次元問題
FEMD 30 プログラム :d.c(6/6) 第一種境界条件 @=0 /* // +---------------------+ // BOUDARY condtons // +---------------------+ */ /* X=Xmn */ =0; js= Inde[]; AMat[jS]= 0.0; Dag[ ]=.0; Rhs [ ]= 0.0; T 0 =0 対角成分 = 右辺 =0 非対角成分 =0 ゼロクリア for(k=0;k<plu;k++){ f(item[k]==0){amat[k]=0.0; 復習 : 一次元問題
FEMD 3 プログラム :d.c(6/6) 第一種境界条件 @=0 /* // +---------------------+ // BOUDARY condtons // +---------------------+ */ /* X=Xmn */ =0; js= Inde[]; AMat[jS]= 0.0; Dag[ ]=.0; Rhs [ ]= 0.0; T 0 =0 対角成分 = 右辺 =0 非対角成分 =0 消去 ゼロクリア for(k=0;k<plu;k++){ f(item[k]==0){amat[k]=0.0; 行列の対称性を保つため 第一種境界条件を適用している節点に対応する 列 を 右辺に移項して消去する ( 今の場合は非対角成分を 0 にするだけで良い ) 復習 : 一次元問題
FEMD 32 第一種境界条件が T 0 の場合 /* // +---------------------+ // BOUDARY condtons // +---------------------+ */ 行列の対称性を保つため 第一種境界条件を適用している節点に対応する 列 を 右辺に移項して消去する /* X=Xmn */ =0; js= Inde[]; AMat[jS]= 0.0; Dag[ ]=.0; Rhs [ ]= PHImn; Dag j j Inde[ j ] k Inde[ j] Amat k Item[ k ] Rhs j for(j=;<;++){ for (k=inde[j];k<inde[j+];k++){ f(item[k]==0){ Rhs [j]= Rhs[j] Amat[k]*PHImn; AMat[k]= 0.0; 復習 : 一次元問題
FEMD 33 第一種境界条件が T 0 の場合 /* // +---------------------+ // BOUDARY condtons // +---------------------+ */ 行列の対称性を保つため 第一種境界条件を適用している節点に対応する 列 を 右辺に移項して消去する /* X=Xmn */ =0; js= Inde[]; AMat[jS]= 0.0; Dag[ ]=.0; Rhs [ ]= PHImn; Dag for(j=;<;++){ for (k=inde[j];k<inde[j+];k++){ f(item[k]==0){ Rhs [j]= Rhs[j] Amat[k]*PHImn; AMat[k]= 0.0; j Rhs Rhs j j j Inde[ j ] k Inde[ j] k k Amat Amat k k s s Amat Item[ k mn s s ] k Item[ k ] where Item[ k s ] 0 復習 : 一次元問題
FEM3D 34 境界条件 :MAT_ASS_BC(2/2) for(n=0;n<;n++){ f( IWKX[n][0] == ){ B[n]= 0.e0; D[n]=.e0; for(k=ndelu[n];k<ndelu[n+];k++){ AMAT[k]= 0.e0; for(n=0;n<;n++){ for(k=ndelu[n];k<ndelu[n+];k++){ f (IWKX[temLU[k]][0] == ) { AMAT[k]= 0.e0; IWKX[n-][0]= となる節点に対して対角成分 = 右辺 =0 非零対角成分 =0 ゼロクリア ここでやっていることも一次元の時と全く変わらない
FEM3D 35 境界条件 :MAT_ASS_BC(2/2) for(n=0;n<;n++){ f( IWKX[n][0] == ){ B[n]= 0.e0; D[n]=.e0; for(k=ndelu[n];k<ndelu[n+];k++){ AMAT[k]= 0.e0; IWKX[n-][0]= となる節点を非零非対角成分として有している節点に対して 右辺へ移項 当該非零非対角成分 =0 for(n=0;n<;n++){ for(k=ndelu[n];k<ndelu[n+];k++){ f (IWKX[temLU[k]][0] == ) { AMAT[k]= 0.e0; 消去 ゼロクリア ここでやっていることも一次元の時と全く変わらない
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 **/ #nclude <stdo.h> #nclude <stdlb.h> FILE* fp_log; #defne GLOBAL_VALUE_DEFIE #nclude "pfem_utl.h" //#nclude "solver.h" etern vod IPUT_CTL(); etern vod IPUT_GRID(); etern vod MAT_CO0(); etern vod MAT_CO(); etern vod MAT_ASS_MAI(); etern vod MAT_ASS_BC(); etern vod SOLVE(); etern vod OUTPUT_UCD(); nt man() { IPUT_CTL(); IPUT_GRID(); MAT_CO0(); MAT_CO(); MAT_ASS_MAI(); MAT_ASS_BC() ; SOLVE(); OUTPUT_UCD() ;