Microsoft PowerPoint - 03-FEM3D-F.ppt [互換モード]

Transcription

1 有限要素法による 三次元定常熱伝導解析プログラム Fortran 編 中島研吾東京大学情報基盤センター

2 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

3 FEM3D 3 対象とする問題 : 三次元定常熱伝導 T T z T z Q z 0 原点から遠い部分が高温 体積当たり発熱量は位置 ( メッシュの中心の座標 ) に依存 Q z C C move

4 FEM3D 有限要素法の処理 支配方程式 ガラーキン法 : 弱形式 要素単位の積分 要素マトリクス生成 全体マトリクス生成 境界条件適用 連立一次方程式

5 FEM3D 5 有限要素法の処理 : プログラム 初期化 制御変数読み込み 座標読み込み 要素生成 (: 節点数 ICELTOT: 要素数 ) 配列初期化 ( 全体マトリクス 要素マトリクス ) 要素 全体マトリクスマッピング (IndeItem) マトリクス生成 要素単位の処理 (do cel= ICELTOT) 要素マトリクス計算 全体マトリクスへの重ね合わせ 境界条件の処理 連立一次方程式 共役勾配法 (CG)

6 FEM3D 6 三次元要素の定式化 三次元熱伝導方程式 ガラーキン法 要素マトリクス生成 プログラムの実行 データ構造 プログラムの構成

7 FEM3D 7 二次元への拡張 : 三角形要素 任意の形状を扱うことができる 特に一次要素は精度が悪く 一部の問題を除いてあまり使用されない

8 FEM3D 二次元への拡張 : 四角形要素 一次元要素と同じ形状関数を 軸に適用することによって 四角形要素の定式化は可能である 三角形と比較して特に低次要素の精度はよい しかしながら 各辺が座標軸に平行な長方形でなければならない 差分法と変わらない 3 このような形状を扱うことができない 2

9 FEM3D 9 アイソパラメトリック要素 (/3) 各要素を 自然座標系 () の正方形要素 [± ±] に変換する 各要素の全体座標系 (global coordnate)() における座標成分を 自然座標系における形状関数 [] ( 従属変数の内挿に使うのと同じ []) を使用して変換する場合 このような要素をアイソパラメトリック要素 (soparametrc element) という

10 アイソパラメトリック要素 (2/3) 2 3 各節点の座標 :( ) ( 2 2 ) ( 3 3 ) ( ) 各節点における温度 :T T 2 T 3 T T T ) ( ) ( ) ( FEM3D 0

11 FEM3D アイソパラメトリック要素 (3/3) 高次の補間関数を使えば 曲線 曲面も扱うことが可能となる そういう意味で 自然座標系 と呼んでいる Sub-Parametrc Super-Parametrc

12 2D 自然座標系の形状関数 (/3) 自然座標系における正方形上の内挿多項式は下式で与えられる : T T T T T T T T T T T T T T T T 3 2 T 各節点での条件より : FEM3D 2

13 2D 自然座標系の形状関数 (2/3) FEM3D T T T T T T T T T T T T T T T T T T T T T T T T T 2 3

14 2D 自然座標系の形状関数 (3/3) 元の式に代入して T について整理すると以下のようになる : ) ( ) ( ) ( ) ( 3 2 形状関数 は以下のようになる : 双一次 (b-lnear) 要素とも呼ばれる 各節点における の値を計算してみよ FEM3D T T T T T

15 FEM3D 5 三次元への拡張 四面体要素 : 二次元における三角形要素 任意の形状を扱うことができる 特に一次要素は精度が悪く 一部の問題を除いてあまり使用されない 高次の四面体要素は広く使用されている 本講義では低次六面体要素 ( アイソパラメトリック要素 ) を使用する tr-lnear 自由度 : 各節点上

16 3D 自然座標系の形状関数 ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( z z T T ) ( ) ( ) ( ) ( FEM3D 6

17 FEM3D 7 三次元要素の定式化 三次元熱伝導方程式 ガラーキン法 要素マトリクス生成 プログラムの実行 データ構造 プログラムの構成

18 ガラーキン法の適用 (/3) 以下のような三次元熱伝導方程式を考慮する ( 熱伝導率一定 ): } { T 要素内の温度分布 ( マトリクス形式 ) 節点における温度を としてある ガラーキン法に従い 重み関数を [] とすると 各要素において以下の積分方程式が得られる : dv Q z T T T T V Q z T T T FEM3D

19 ガラーキン法の適用 (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 dv T T T dv T T T V z T z T T zz T V FEM3D 9

20 FEM3D 20 ガラーキン法の適用 (3/3) 体積あたり発熱量の項を加えて次式が得られる : Q T T T z z dv Q dv 0 V この式を弱形式 (weak form) と呼ぶ 元の微分方程式では 2 階の微分が含まれていたが 上式では グリーンの定理によって 階微分に低減されている 弱形式によって近似関数 ( 形状関数 内挿関数 ) に対する要求が弱くなっている : すなわち線形関数で2 階微分の効果を記述できる 項が増えただけで 一次元と同じ V

21 境界条件を考慮した弱形式 : 各要素 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 2

22 FEM3D 22 要素マトリクス : 行列 j j k j

23 FEM3D 23 要素マトリクス :k j j k j dv T ( e) k T V V T z z dv V dv k j V j j z j z dv

24 D-Part2 2 数値的に積分を実施する方法 台形公式 シンプソンの公式 ガウスの積分公式 ( ガウス=ルジャンドル (Gauss-Legendre) とも呼ばれる 精度良い ) 不定積分を解析的に求めるのではなく 有限な数のサンプル点における値を利用する X X 2 f m d wk f k k

25 D-Part2 25 ガウス積分公式 : 一次元の例シンプソンの公式より精度良い Smpson s f Gauss ( ) sn( ) X h S X X 2 X 3 0 X 2 X 3 2 X h 3 2 X X 3 X f X f X f X S h 2 2 / k f d W k f k f h d

26 D-Part2 26 ガウスの積分公式 無次元化された自然座標系 [-+] を対象とする m 個の積分点を使用すると (2m-) 次の関数までは近似可能 ( 従って 次 2 次の内挿関数 ( 形状関数 ) を使用するときは m=2で十分) f m m 2 m 3 m d w k f k k k 0.00 w 2.00 k wk.00 k 0.00 wk / w 5/ 9 k k k =- =0 =+ k = k =

27 Gaussan Quadrature ガウスの積分公式 自然座標系 () 上で定義 ガウス積分公式が使える ( 三次元 ) 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 27

28 2 Gaussan Quadrature ガウスの積分公式 この組み合わせがよく 使われる 二次元では 点におけるfの値を計算 して足せば良い 三次 元では点

29 29 Gaussan Quadrature ガウスの積分公式 この組み合わせがよく 使われる 二次元では 点におけるfの値を計算 して足せば良い 三次 元では点 I f ( ) d d W W f ( ) m n j j j.0.0 f ( ).0.0 f ( ).0.0 f ( ).0.0 f ( )

30 あとは積分を求めれば良い 自然座標系 () 上で定義 ガウス積分公式が使える ( 三次元 ) しかし 微分が 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 30

31 自然座標系における偏微分 (/) 偏微分の公式より以下のようになる : z z z z z z ) ( ) ( ) ( z は定義より簡単に求められるがを実際の計算で使用する FEM3D 3

32 自然座標系における偏微分 (2/) マトリックス表示すると : z J z z z z J : ヤコビのマトリクス (Jacob matr Jacoban) FEM3D 32

33 自然座標系における偏微分 (3/) の定義より簡単に求められる z z z J J J z z z J J J z z z J J J FEM3D 33

34 自然座標系における偏微分 (/) 従って下記のように偏微分を計算できる ヤコビアン (3 3 行列 ) の逆行列を求める J z z z z FEM3D 3

35 要素単位での積分 V j j j V z j z j j j dv z z dv k j k j FEM3D 35

36 自然座標系での積分 d d d J z z dddz z z dv z z j j j j j j V j j j det j k j FEM3D 36

37 ガウスの積分公式 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 37 j k j

38 FEM3D 3 残りの手順 ここまでで 要素ごとの積分が可能となる あとは : 全体マトリクスへの重ね合わせ 境界条件処理 連立一次方程式を解く 詳細はプログラムの内容を解説しながら説明する

39 要素 全体マトリクス重ね合わせ () () 3 () 2 () () () 3 () 2 () () () 3 () 2 () () 3 () 33 () 32 () 3 () 2 () 23 () 22 () 2 () () 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) (2) 3 (2) 2 (2) (2) (2) 3 (2) 2 (2) (2) (2) 3 (2) 2 (2) (2) 3 (2) 33 (2) 32 (2) 3 (2) 2 (2) 23 (2) 22 (2) 2 (2) (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 } { } ]{ [ 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 39

40 FEM3D 0 三次元要素の定式化 三次元弾性力学方程式 ガラーキン法 要素マトリクス生成 演習 プログラムの実行 データ構造 プログラムの構成

41 FEM3D 演習 ガウスの積分公式を使用して以下の四角形の面積を求めよ ( プログラムを作って 計算機で計算してください ) V 3 2 :(.0.0) 2:(.0 2.0) 3:( ) :(2.0.0) I V dv det J dd

42 FEM3D 2 ヒント (/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()= d0 POI(2)= d0 SUM= 0.d0 do jp= 2 do p= 2 FC = F(POI(p)POI(jp)) SUM= SUM + W (p)*w (jp)*fc

43 ヒント (2/2) J J det ) ( ) ( ) ( ) ( 3 2 FEM3D 3

44 FEM3D 三次元要素の定式化 三次元弾性力学方程式 ガラーキン法 要素マトリクス生成 プログラムの実行 データ構造 プログラムの構成

45 FEM3D 5 対象とする問題 : 三次元定常熱伝導 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

46 FEM3D 6 インストール (Cgwn では *.ee) インストール >$ cd <$P-TOP>/fem3d/src >$ make >$ ls../run/sol sol メッシュジェネレータインストール >$ cd <$P-TOP>/fem3d/run >$ gfortran O3 mgcube.f o mgcube

47 FEM3D 7 計算の流れメッシュ生成 計算 ファイル名称固定 mgcube メッシュジェネレータ cube.0 メッシュファイル sol FEM 本体 IPUT.DAT 制御データ test.np 可視化用出力

48 FEM3D メッシュ生成 (Cgwn では mgcube.ee) Z >$ cd <$P-TOP>/fem3d/run >$./mgcube T=0@Z=z ma X Y Z 各辺長さを訊いてくる このように入れてみる Z >$ ls cube.0 生成を確認 cube.0 X Y X Y

49 FEM3D 9 制御ファイル :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

50 FEM3D 50 実行 (Cgwn では sol.ee >$ cd <$P-TOP>/fem3d/run >$./sol >$ ls test.np 生成を確認 test.np

51 FEM3D 5 ParaVew ファイルを開く 図の表示 イメージファイルの保存

52 FEM3D 52 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

53 FEM3D 53 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.

54 FEM3D 5 UCD Format (3/3): Old Format ( 全節点数 ) ( 全要素数 ) ( 各節点のデータ数 ) ( 各要素のデータ数 ) ( モデルのデータ数 ) ( 節点番号 ) (X 座標 ) (Y 座標 ) (Z 座標 ) ( 節点番号 2) (X 座標 ) (Y 座標 ) (Z 座標 ) ( 要素番号 ) ( 材料番号 ) ( 要素の種類 ) ( 要素を構成する節点のつながり ) ( 要素番号 2) ( 材料番号 ) ( 要素の種類 ) ( 要素を構成する節点のつながり ) ( 要素のデータ成分数 ) ( 成分 の構成数 ) ( 成分 2 の構成数 ) ( 各成分の構成数 ) ( 要素データ成分 のラベル )( 単位 ) ( 要素データ成分 2 のラベル )( 単位 ) ( 各要素データ成分のラベル )( 単位 ) ( 要素番号 ) ( 要素データ ) ( 要素データ 2) ( 要素番号 2) ( 要素データ ) ( 要素データ 2) ( 節点のデータ成分数 ) ( 成分 の構成数 ) ( 成分 2 の構成数 ) ( 各成分の構成数 ) ( 節点データ成分 のラベル )( 単位 ) ( 節点データ成分 2 のラベル )( 単位 ) ( 各節点データ成分のラベル )( 単位 ) ( 節点番号 ) ( 節点データ ) ( 節点データ 2) ( 節点番号 2) ( 節点データ ) ( 節点データ 2)

55 FEM3D 55 三次元要素の定式化 三次元弾性力学方程式 ガラーキン法 要素マトリクス生成 プログラムの実行 データ構造 プログラムの構成

56 FEM3D 56 メッシュファイル構成 :cube.0 番号は から始まっている 節点データ 節点数 節点番号 座標 要素データ 要素数 要素タイプ 要素番号 材料番号 コネクティビティ 節点グループデータ グループ数 グループ内節点数 グループ名 グループ内節点

57 FEM3D 57 cube.0: 節点データ (X=Y=Z=) 25 =5*5*5( 節点数 ) Z T=0@Z=z ma Z X Y X Y 節点 ID X 座標 Y 座標 Z 座標 move

58 FEM3D 5 cube.0: 要素データ (/2) 6 =**( 要素数 ) Z T=0@Z=z ma Z 要素タイプ :36 三次元 六面体 一次 Y X Y X move

59 FEM3D 59 cube.0: 要素データ (2/2) 要素 ID 材料 ID 節点 ID X Z T=0@Z=z ma Z Y X Y

60 FEM3D 60 X=Y=Z= XP=YP=ZP=5 ICELTOT= 6 IODTOT= 25 IBODTOT= 25 k= k= k= k= k=

61 FEM3D 6 cube.0: 節点グループデータ グループ総数 各グループ節点数 ( 累積 ) Z Xmn T=0@Z=z ma Ymn Z Zmn Zma X Y X Y ( 以下使用せず )

62 FEM3D 62 X=Y=Z= XP=YP=ZP=5 ICELTOT= 6 IODTOT= 25 IBODTOT= 25 k= k= k= k= k= Xmn: = Ymn: j= Zmn: k= Zma:k=5

63 FEM3D 63 メッシュ生成 実は技術的には大きな課題 複雑形状 大規模メッシュ 並列化が難しい 市販のメッシュ生成アプリケーション FEMAP CAD データとのインタフェース move

64 FEM3D 6 三次元要素の定式化 三次元弾性力学方程式 ガラーキン法 要素マトリクス生成 プログラムの実行 データ構造 プログラムの構成

65 FEM3D 65 有限要素法の処理 : プログラム 初期化 制御変数読み込み 座標読み込み 要素生成 (: 節点数 E: 要素数 ) 配列初期化 ( 全体マトリクス 要素マトリクス ) 要素 全体マトリクスマッピング (IndeItem) マトリクス生成 要素単位の処理 (do cel= E) 要素マトリクス計算 全体マトリクスへの重ね合わせ 境界条件の処理 連立一次方程式 共役勾配法 (CG)

66 FEM3D 66 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 法計算

67 FEM3D 67 全体処理 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

68 FEM3D 6 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 右辺ベクトル 未知数ベクトル

69 FEM3D 69 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 諸定数 ( 実数 )

70 FEM3D 70 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 熱伝導率 体積当たり発熱量係数

71 FEM3D 7 制御ファイル入力 :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

72 FEM3D 72 メッシュ入力 :IPUT_GRID(/3)!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)

73 FEM3D 73 cube.0: 節点データ (X=Y=Z=) 25 = Z T=0@Z=z ma Z X Y X Y XYZ(3) move

74 FEM3D 7 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 並みに簡単にやるための関数

75 FEM3D 75 メッシュ入力 :IPUT_GRID(2/3)!C!C-- ELEMET read (*) ICELTOT allocate (ICELOD(ICELTOT)) read ('(00)') (TYPE = ICELTOT) do cel= ICELTOT read ('(0025)') IMAT (ICELOD(celk) k=)

76 FEM3D 76 cube.0: 要素データ (/2) 6 = ICELTOT Z T=0@Z=z ma Z 要素タイプ :36 三次元 六面体 一次 Y X Y X move

77 FEM3D 77 cube.0: 要素データ (2/2) MAT ICELOD(cel) X Z T=0@Z=z ma Z Y X Y

78 FEM3D 7 メッシュ入力 :IPUT_GRID(3/3)!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

79 FEM3D 79 cube.0: 節点グループデータ ODGRPtot ODGTP_IDEX(-) Xmn ODGRP_AME() ODGRP_ITEM(-25) Ymn ODGRP_AME(2) ODGRP_ITEM(26-50) Zmn ODGRP_AME(3) ODGRP_ITEM(5-75) Zma ODGRP_AME() ODGRP_ITEM(76-00)

80 FEM3D 0 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 法計算

81 FEM3D 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 右辺ベクトル 未知数ベクトル

82 FEM3D 2 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 諸定数 ( 実数 )

83 FEM3D 3 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 熱伝導率 体積当たり発熱量係数

84 FEM3D マトリクス生成まで 一次元のときは ndetem に関連した情報を簡単に作ることができた 非ゼロ非対角成分の数は2 番号が自分に対して :+と- 三次元の場合はもっと複雑 非ゼロ非対角成分の数は7~26( 現在の形状 ) 実際はもっと複雑 前以て 非ゼロ非対角成分数はわからない move

85 FEM3D 5 マトリクス生成まで 一次元のときは ndetem に関連した情報を簡単に作ることができた 非ゼロ非対角成分の数は2 番号が自分に対して :+と- 三次元の場合はもっと複雑 非ゼロ非対角成分の数は7~26( 現在の形状 ) 実際はもっと複雑 前以て 非ゼロ非対角成分の数はわからない ILU()IALU(LU) を使って非ゼロ非対角成分数を予備的に勘定する

86 FEM3D 6 全体処理 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 MAT_CO0: IU IALU 生成 MAT_CO: nde tem 生成 とりあえず から始まる節点番号を記憶

87 FEM3D 7 MAT_CO0: 全体構成 do cel= ICELTOT 節点相互の関係から ILUIALU を生成 (FID_ODE)

88 FEM3D 行列コネクティビティ生成 : MAT_CO0(/)!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: 各節点における非ゼロ非対角成分の最大数 ( 接続する節点数 ) 今の問題の場合はわかっているので このようにできる 不明の場合の実装 : レポート課題

89 FEM3D 9 行列コネクティビティ生成 : MAT_CO0(/)!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) 各節点の非零非対角成分数 各節点の非零非対角成分 ( 列番号 )

90 FEM3D 90 行列コネクティビティ生成 : MAT_CO0(2/): から始まる番号 do cel= ICELTOT n= ICELOD(cel) n2= ICELOD(cel2) n3= ICELOD(cel3) n= ICELOD(cel) 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 (nn) call FID_TS_ODE (nn5) call FID_TS_ODE (nn6) call FID_TS_ODE (nn7) call FID_TS_ODE (nn) call FID_TS_ODE (n2n) call FID_TS_ODE (n2n3) call FID_TS_ODE (n2n) 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 (n3n) call FID_TS_ODE (n3n5) call FID_TS_ODE (n3n6) call FID_TS_ODE (n3n7) call FID_TS_ODE (n3n)

91 FEM3D 9 節点探索 :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) 各節点の非零非対角成分数 各節点の非零非対角成分 ( 列番号 )

92 FEM3D 92 節点探索 :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

93 FEM3D 93 節点探索 :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 に格納

94 FEM3D 9 行列コネクティビティ生成 : MAT_CO0(3/) call FID_TS_ODE (nn) call FID_TS_ODE (nn2) call FID_TS_ODE (nn3) call FID_TS_ODE (nn5) call FID_TS_ODE (nn6) call FID_TS_ODE (nn7) call FID_TS_ODE (nn) call FID_TS_ODE (n5n) call FID_TS_ODE (n5n2) call FID_TS_ODE (n5n3) call FID_TS_ODE (n5n) 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 (n6n) call FID_TS_ODE (n6n5) call FID_TS_ODE (n6n7) call FID_TS_ODE (n6n) call FID_TS_ODE (n7n) call FID_TS_ODE (n7n2) call FID_TS_ODE (n7n3) call FID_TS_ODE (n7n) call FID_TS_ODE (n7n5) call FID_TS_ODE (n7n6) call FID_TS_ODE (n7n)

95 FEM3D 95 行列コネクティビティ生成 : MAT_CO0(/) call FID_TS_ODE (nn) call FID_TS_ODE (nn2) call FID_TS_ODE (nn3) call FID_TS_ODE (nn) 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 程度のものをソートする

96 FEM3D 96 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

97 FEM3D 97 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

98 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) tem に から始まる節点番号を記憶 deallocate (ILU IALU) end subroutne MAT_CO

99 FEM3D 99 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 これらはもはや不要

100 FEM3D 00 全体処理 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

101 FEM3D 0 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

102 FEM3D 02 係数行列 :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()= D+00 WEI(2)= D+00 POS()= D+00 POS(2)= D+00

103 FEM3D 03 係数行列 :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()= D+00 WEI(2)= D+00 POS()= D+00 POS(2)= D+00 POS: WEI: 積分点座標重み係数

104 FEM3D 0 係数行列 :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(pjpkp)= Oth * QM * EP * TM SHAPE(pjpkp5)= Oth * QM * EM * TP SHAPE(pjpkp6)= Oth * QP * EM * TP SHAPE(pjpkp7)= Oth * QP * EP * TP SHAPE(pjpkp)= Oth * QM * EP * TP

105 FEM3D 05 係数行列 :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(pjpkp)= Oth * QM * EP * TM SHAPE(pjpkp5)= Oth * QM * EM * TP SHAPE(pjpkp6)= Oth * QP * EM * TP SHAPE(pjpkp7)= Oth * QP * EP * TP SHAPE(pjpkp)= Oth * QM * EP * TP QP EP TP k QM j EM j j TMk k k

106 FEM3D 06 係数行列 :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(pjpkp)= Oth * QM * EP * TM SHAPE(pjpkp5)= Oth * QM * EM * TP SHAPE(pjpkp6)= Oth * QP * EM * TP SHAPE(pjpkp7)= Oth * QP * EP * TP SHAPE(pjpkp)= Oth * QM * EP * TP

107 FEM3D 07 係数行列 :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(pjpkp)= Oth * QM * EP * TM SHAPE(pjpkp5)= Oth * QM * EM * TP SHAPE(pjpkp6)= Oth * QP * EM * TP SHAPE(pjpkp7)= Oth * QP * EP * TP SHAPE(pjpkp)= Oth * QM * EP * TP ) ( ) ( ) ( ) ( 3 2 ) ( ) ( ) ( ) ( 7 6 5

108 FEM3D 0 係数行列 :MAT_ASS_MAI(3/6) PQ(jpkp)= - Oth * EM * TM PQ(jpkp2)= + Oth * EM * TM PQ(jpkp3)= + Oth * EP * TM PQ(jpkp)= - 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(pkp)= + 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(pjp)= - 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) n= ICELOD(cel) n5= ICELOD(cel5) n6= ICELOD(cel6) n7= ICELOD(cel7) n= ICELOD(cel) ( j k ) l PQ( j k) ( j k ) PE( k) l ( j k ) PT ( l k j) ( j ) ( j k ) 2 ( j k ) 3 ( j k ) 3 ( j k ) j j j j における形状関数の一階微分 k k k k

109 FEM3D 09 係数行列 :MAT_ASS_MAI(3/6) PQ(jpkp)= - Oth * EM * TM PQ(jpkp2)= + Oth * EM * TM PQ(jpkp3)= + Oth * EP * TM PQ(jpkp)= - 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(pkp)= + 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(pjp)= - 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) n= ICELOD(cel) n5= ICELOD(cel5) n6= ICELOD(cel6) n7= ICELOD(cel7) n= ICELOD(cel)

110 FEM3D 0 係数行列 :MAT_ASS_MAI(/6) nodlocal()= n nodlocal(2)= n2 nodlocal(3)= n3 nodlocal()= n nodlocal(5)= n5 nodlocal(6)= n6 nodlocal(7)= n7 nodlocal()= n 節点の節点番号 X= XYZ(n) X2= XYZ(n2) X3= XYZ(n3) X= XYZ(n) X5= XYZ(n5) X6= XYZ(n6) X7= XYZ(n7) X= XYZ(n) Y= XYZ(n2) Y2= XYZ(n22) Y3= XYZ(n32) Y= XYZ(n2) Y5= XYZ(n52) Y6= XYZ(n62) Y7= XYZ(n72) Y= XYZ(n2) QVC= Oth * (X+X2+X3+X+X5+X6+X7+X+ & Y+Y2+Y3+Y+Y5+Y6+Y7+Y) Z= XYZ(n3) Z2= XYZ(n23) Z3= XYZ(n33) Z= XYZ(n3) Z5= XYZ(n53) Z6= XYZ(n63) Z7= XYZ(n73) Z= XYZ(n3) call JACOBI (DETJ PQ PE PT PX PY PZ & & X X2 X3 X X5 X6 X7 X & & Y Y2 Y3 Y Y5 Y6 Y7 Y & & Z Z2 Z3 Z Z5 Z6 Z7 Z )

111 FEM3D 係数行列 :MAT_ASS_MAI(/6) nodlocal()= n nodlocal(2)= n2 nodlocal(3)= n3 nodlocal()= n nodlocal(5)= n5 nodlocal(6)= n6 nodlocal(7)= n7 nodlocal()= n X= XYZ(n) X2= XYZ(n2) X3= XYZ(n3) X= XYZ(n) X5= XYZ(n5) X6= XYZ(n6) X7= XYZ(n7) X= XYZ(n) Y= XYZ(n2) Y2= XYZ(n22) Y3= XYZ(n32) Y= XYZ(n2) Y5= XYZ(n52) Y6= XYZ(n62) Y7= XYZ(n72) Y= XYZ(n2) QVC= Oth * (X+X2+X3+X+X5+X6+X7+X+ & Y+Y2+Y3+Y+Y5+Y6+Y7+Y) Z= XYZ(n3) Z2= XYZ(n23) Z3= XYZ(n33) Z= XYZ(n3) Z5= XYZ(n53) Z6= XYZ(n63) Z7= XYZ(n73) Z= XYZ(n3) 節点の X 座標 節点の Y 座標 節点の Z 座標 call JACOBI (DETJ PQ PE PT PX PY PZ & & X X2 X3 X X5 X6 X7 X & & Y Y2 Y3 Y Y5 Y6 Y7 Y & & Z Z2 Z3 Z Z5 Z6 Z7 Z )

112 FEM3D 2 係数行列 :MAT_ASS_MAI(/6) nodlocal()= n nodlocal(2)= n2 nodlocal(3)= n3 nodlocal()= n nodlocal(5)= n5 nodlocal(6)= n6 nodlocal(7)= n7 nodlocal()= n X= XYZ(n) X2= XYZ(n2) X3= XYZ(n3) X= XYZ(n) X5= XYZ(n5) X6= XYZ(n6) X7= XYZ(n7) X= XYZ(n) Y= XYZ(n2) Y2= XYZ(n22) Y3= XYZ(n32) Y= XYZ(n2) Y5= XYZ(n52) Y6= XYZ(n62) Y7= XYZ(n72) Y= XYZ(n2) QVC= Oth * (X+X2+X3+X+X5+X6+X7+X+ & Y+Y2+Y3+Y+Y5+Y6+Y7+Y) Z= XYZ(n3) Z2= XYZ(n23) Z3= XYZ(n33) Z= XYZ(n3) Z5= XYZ(n53) Z6= XYZ(n63) Z7= XYZ(n73) Z= XYZ(n3) 節点の X 座標 節点の Y 座標 T Q T 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 X X5 X6 X7 X & & Y Y2 Y3 Y Y5 Y6 Y7 Y & & Z Z2 Z3 Z Z5 Z6 Z7 Z )

113 FEM3D 3 係数行列 :MAT_ASS_MAI(/6) nodlocal()= n nodlocal(2)= n2 nodlocal(3)= n3 nodlocal()= n nodlocal(5)= n5 nodlocal(6)= n6 nodlocal(7)= n7 nodlocal()= n X= XYZ(n) X2= XYZ(n2) X3= XYZ(n3) X= XYZ(n) X5= XYZ(n5) X6= XYZ(n6) X7= XYZ(n7) X= XYZ(n) Y= XYZ(n2) Y2= XYZ(n22) Y3= XYZ(n32) Y= XYZ(n2) Y5= XYZ(n52) Y6= XYZ(n62) Y7= XYZ(n72) Y= XYZ(n2) QVC= Oth * (X+X2+X3+X+X5+X6+X7+X+ & Y+Y2+Y3+Y+Y5+Y6+Y7+Y) Z= XYZ(n3) Z2= XYZ(n23) Z3= XYZ(n33) Z= XYZ(n3) Z5= XYZ(n53) Z6= XYZ(n63) Z7= XYZ(n73) Z= XYZ(n3) T Q T 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 X X5 X6 X7 X & & Y Y2 Y3 Y Y5 Y6 Y7 Y & & Z Z2 Z3 Z Z5 Z6 Z7 Z )

114 FEM3D 係数行列 :MAT_ASS_MAI(/6) nodlocal()= n nodlocal(2)= n2 nodlocal(3)= n3 nodlocal()= n nodlocal(5)= n5 nodlocal(6)= n6 nodlocal(7)= n7 nodlocal()= n X= XYZ(n) X2= XYZ(n2) X3= XYZ(n3) X= XYZ(n) X5= XYZ(n5) X6= XYZ(n6) X7= XYZ(n7) X= XYZ(n) Y= XYZ(n2) Y2= XYZ(n22) Y3= XYZ(n32) Y= XYZ(n2) Y5= XYZ(n52) Y6= XYZ(n62) Y7= XYZ(n72) Y= XYZ(n2) QVC= Oth * (X+X2+X3+X+X5+X6+X7+X+ & Y+Y2+Y3+Y+Y5+Y6+Y7+Y) Z= XYZ(n3) Z2= XYZ(n23) Z3= XYZ(n33) Z= XYZ(n3) Z5= XYZ(n53) Z6= XYZ(n63) Z7= XYZ(n73) Z= XYZ(n3) call JACOBI (DETJ PQ PE PT PX PY PZ & & X X2 X3 X X5 X6 X7 X & & Y Y2 Y3 Y Y5 Y6 Y7 Y & & Z Z2 Z3 Z Z5 Z6 Z7 Z )

115 FEM3D 5 JACOBI(/) subroutne JACOBI (DETJ PQ PE PT PX PY PZ & & X X2 X3 X X5 X6 X7 X Y Y2 Y3 Y Y5 Y6 Y7 Y & & Z Z2 Z3 Z 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(pjpkp)=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(pjpkp)=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(pjpkp)=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

116 FEM3D 6 自然座標系における偏微分 (/) 偏微分の公式より以下のようになる : z z z z z z ) ( ) ( ) ( z は定義より簡単に求められるがを実際の計算で使用する

117 FEM3D 7 自然座標系における偏微分 (2/) マトリックス表示すると : z J z z z z J : ヤコビのマトリクス (Jacob matr Jacoban) J J J J J J J J J J

118 FEM3D 自然座標系における偏微分 (3/) の定義より簡単に求められる z z z J J J z z z J J J z z z J J J

119 FEM3D 9!C!C== DETERMIAT of the JACOBIA JACOBI(2/) dxdq = & & + PQ(jpkp) * X + PQ(jpkp2) * X2 & & + PQ(jpkp3) * X3 + PQ(jpkp) * X & & + PQ(jpkp5) * X5 + PQ(jpkp6) * X6 & & + PQ(jpkp7) * X7 + PQ(jpkp) * X dydq = & & + PQ(jpkp) * Y + PQ(jpkp2) * Y2 & & + PQ(jpkp3) * Y3 + PQ(jpkp) * Y & 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(jpkp) * Z & & + 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(pkp) * X & & + 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(pkp) * Y & & + 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(pkp) * Z & & + PE(pkp5) * Z5 + PE(pkp6) * Z6 & z & + PE(pkp7) * Z7 + PE(pkp) * Z dzdq J3 J J J

120 FEM3D 20 JACOBI(3/) dxdt = & & + PT(pjp) * X + PT(pjp2) * X2 & & + PT(pjp3) * X3 + PT(pjp) * X & & + PT(pjp5) * X5 + PT(pjp6) * X6 & & + PT(pjp7) * X7 + PT(pjp) * X dydt = & & + PT(pjp) * Y + PT(pjp2) * Y2 & & + PT(pjp3) * Y3 + PT(pjp) * Y & & + PT(pjp5) * Y5 + PT(pjp6) * Y6 & & + PT(pjp7) * Y7 + PT(pjp) * Y dzdt = & & + PT(pjp) * Z + PT(pjp2) * Z2 & & + PT(pjp3) * Z3 + PT(pjp) * Z & & + 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 J J J a3= coef * ( dxde*dydt - dyde*dxdt ) a32= coef * ( dydq*dxdt - dxdq*dydt ) a33= coef * ( dxdq*dyde - dydq*dxde ) DETJ(pjpkp)= dabs(detj(pjpkp))

121 FEM3D 2 自然座標系における偏微分 (/) 従って下記のように偏微分を計算できる ヤコビアン (3 3 行列 ) の逆行列を求める J z z z z

122 FEM3D 22 JACOBI(3/) dxdt = & & + PT(pjp) * X + PT(pjp2) * X2 & & + PT(pjp3) * X3 + PT(pjp) * X & & + PT(pjp5) * X5 + PT(pjp6) * X6 & & + PT(pjp7) * X7 + PT(pjp) * X dydt = & & + PT(pjp) * Y + PT(pjp2) * Y2 & & + PT(pjp3) * Y3 + PT(pjp) * Y & & + PT(pjp5) * Y5 + PT(pjp6) * Y6 & & + PT(pjp7) * Y7 + PT(pjp) * Y dzdt = & & + PT(pjp) * Z + PT(pjp2) * Z2 & & + PT(pjp3) * Z3 + PT(pjp) * Z & & + 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 a a a DETJ(pjpkp)= dabs(detj(pjpkp))

123 FEM3D 23 JACOBI(/)!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(pjpkp)= a*pq(jpkp) + a2*pe(pkp) + a3*pt(pjp) 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(pjpkp)= a2*pq(jpkp) + a22*pe(pkp) + a23*pt(pjp) 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(pjpkp)= a3*pq(jpkp) + a32*pe(pkp) + a33*pt(pjp) 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

124 FEM3D 2 係数行列 :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 全体行列の非対角成分 A p jp kk:tem におけるアドレス p= nodlocal(e) jp= nodlocal(je) 2

125 FEM3D 25 要素マトリクス : 行列 j j k j

126 FEM3D 26 係数行列 :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 におけるアドレス

127 FEM3D 27 係数行列 :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

128 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) 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 I L j M j z f ( ) ddd W W j Wk f ( j k ) k z j det J ddd

129 FEM3D 29 係数行列 :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) 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 I L j M j z f ( ) ddd W W j Wk f ( j k ) k z j det J ddd

130 FEM3D 30 係数行列 :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

131 FEM3D 3 係数行列 :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

132 FEM3D 32 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

133 FEM3D 33 境界条件 :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)= とする

134 FEM3D 3 境界条件 :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

135 FEM3D 35 対象とする問題 : 一次元熱伝導問題 体積当たり一様発熱 Q T Q 0 =0 ( mn ) = ma 一様な : 断面積 A 熱伝導率 体積当たり一様発熱 ( 時間当たり ) QL -3 T - 境界条件 =0 := 0 ( 固定 ) T = ma : 0( 断熱 ) 復習 : 一次元問題 Q

136 FEM3D 36 =0 で成立する方程式 T =0 体積当たり一様発熱 Q T Q 0 =0 ( mn ) = ma 一様な : 断面積 A 熱伝導率 体積当たり一様発熱 ( 時間当たり ) QL -3 T - 境界条件 =0 :=0 ( 固定 ) T = ma : 0( 断熱 ) Q 復習 : 一次元問題

137 FEM3D 37 プログラム :d.f(6/6) !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 復習 : 一次元問題

138 FEM3D 3 プログラム :d.f(6/6) !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 復習 : 一次元問題

139 FEM3D 39 プログラム :d.f(6/6) !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 にするだけで良い ) 復習 : 一次元問題

140 FEM3D 0 第一種境界条件が 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 復習 : 一次元問題

141 FEM3D 第一種境界条件が 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 ) 復習 : 一次元問題

142 FEM3D 2 境界条件 :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 ここでやっていることも一次元の時と全く変わらない

143 FEM3D 3 境界条件 :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 消去 ゼロクリア ここでやっていることも一次元の時と全く変わらない

144 FEM3D 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 法計算

145 FEM3D 5 全体処理 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

146 FEM3D 6 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

147 FEM3D 7 前処理付き共役勾配法 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 (-)

148 FEM3D 対角スケーリング 点ヤコビ前処理 前処理行列として もとの行列の対角成分のみを取り出した行列を前処理行列 [M] とする 対角スケーリング 点ヤコビ (pont-jacob) 前処理 M D D D D solve [M]z (-) = r (-) という場合に逆行列を簡単に求めることができる

149 FEM3D 9 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= nteger(knd=knt ) :: MAXIT real (knd=kreal) :: TOL W SS

150 FEM3D 50 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()= 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= 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 (-)

151 FEM3D 5 CG ソルバー (2/6)!C!C !C IIT.!C !C=== ERROR= 0 allocate (WW()) 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()= WW(DD) /{D} 対角成分の逆数 ( 前処理用 ) その都度 除算をすると効率が悪いため 予め配列に格納

152 FEM3D 52 CG ソルバー (2/6)!C!C !C IIT.!C !C=== ERROR= 0 allocate (WW()) 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 (-)

153 FEM3D 53 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===

154 FEM3D 5 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

155 FEM3D 55 CG ソルバー (/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 (-)

156 FEM3D 56 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= C0 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 (-)

157 FEM3D 57 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 (-)

158 FEM3D 5!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