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

Transcription

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

2 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

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

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 自然座標系の形状関数 ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( FEM3D 6 T T ) ( ) ( ) ( ) (

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

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

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

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

21 境界条件を考慮した弱形式 : 各要素 dv dv dv V T V T V T e ) ( ) ( ) ( ) ( e e e f dv Q f V T e ) ( FEM3D 2

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

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

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

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 / f d W f f h d

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

27 Gaussan Quadrature ガウスの積分公式 自然座標系 () 上で定義 ガウス積分公式が使える ( 三次元 ) j W W W ) ( j LM: 方向の積分点数 : 積分点での重み係数 : 積分点の座標値 j 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 あとは積分を求めれば良い 自然座標系 () 上で定義 ガウス積分公式が使える ( 三次元 ) しかし 微分が j W W W ) ( j LM: 方向の積分点数 : 積分点での重み係数 : 積分点の座標値 j j M j L f W W W d d d f I ) ( ) ( FEM3D 30

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

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

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

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

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

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

37 ガウスの積分公式 j j M j L f W W W d d d f I ) ( ) ( d d d J j j j det FEM3D 37 j 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 f (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 f { ]{ [ 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 enddo enddo

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

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

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

46 FEM3D 6 インストール (Cgwn では *.ee) インストール >$ cd <$P-TOP>/fem3d/src >$ mae >$ ls../run/sol sol メッシュジェネレータインストール >$ cd <$P-TOP>/fem3d/run >$ cc O3 mgcube.c 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= 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 QVOL C C T Q 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= ma Z X Y X Y 節点 ID X 座標 Y 座標 Z 座標 move

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

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

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

61 FEM3D 6 cube.0: 節点グループデータ グループ総数 各グループ節点数 ( 累積 ) Z Xmn T=0@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 = = = = = Xmn: = Ymn: j= Zmn: = Zma:=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 **/ #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() ;

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

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

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

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

72 FEM3D 72 メッシュ入力 :IPUT_GRID(/3) #nclude <stdo.h> #nclude <stdlb.h> #nclude "pfem_utl.h" #nclude "allocate.h" vod IPUT_GRID() { FILE *fp; nt jnncelse; 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]);

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

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

75 FEM3D 75 メッシュ入力 :IPUT_GRID(2/3) /** **/ ELEMET fscanf(fp"%d"&iceltot); ICELOD=(KIT**)allocate_matr(seof(KIT)ICELTOT); for(=0;<iceltot;++) fscanf(fp"%d"&type); ICELOD[][j] の中身としては から始まる通し節点番号がそのまま読み込まれている 要素番号は 0 から番号付け 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][]&ICELOD[cel][5]&ICELOD[cel][6]&ICELOD[cel][7]);

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

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

78 FEM3D 7 メッシュ入力 :IPUT_GRID(3/3) /** **/ ODE grp. nfo. fscanf(fp"%d"&odgrptot); 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(=0;<odgrptot;++){ S= ODGRP_IDEX[]; E= ODGRP_IDEX[+]; fscanf(fp"%s"odgrp_ame[].name); nn= E - S; f( nn!= 0 ){ for(=s;<e;++) fscanf(fp"%d"&odgrp_item[]); fclose(fp); 節点グループの中身も から始まる通し節点番号がそのまま読み込まれている

79 FEM3D 79 cube.0: 節点グループデータ ODGRPtot ODGTP_IDEX[-] Xmn ODGRP_AME[] ODGRP_ITEM[0-2] Ymn ODGRP_AME[2] ODGRP_ITEM[25-9] Zmn ODGRP_AME[3] ODGRP_ITEM[50-7] Zma ODGRP_AME[] ODGRP_ITEM[75-99]

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

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

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

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

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

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

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

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

89 FEM3D 9 行列コネクティビティ生成 : MAT_CO0(/) /** ** 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 jceln; nt nn2n3nn5n6n7n; 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;

90 FEM3D 90 行列コネクティビティ生成 : MAT_CO0(2/): から始まる番号 for( cel=0;cel< ICELTOT;cel++){ n=icelod[cel][0]; n2=icelod[cel][]; n3=icelod[cel][2]; n=icelod[cel][3]; n5=icelod[cel][]; n6=icelod[cel][5]; n7=icelod[cel][6]; n=icelod[cel][7]; FID_TS_ODE (nn2); FID_TS_ODE (nn3); FID_TS_ODE (nn); FID_TS_ODE (nn5); FID_TS_ODE (nn6); FID_TS_ODE (nn7); FID_TS_ODE (nn); FID_TS_ODE (n2n); FID_TS_ODE (n2n3); FID_TS_ODE (n2n); 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 (n3n); FID_TS_ODE (n3n5); FID_TS_ODE (n3n6); FID_TS_ODE (n3n7); FID_TS_ODE (n3n);

91 FEM3D 9 節点探索 :FID_TS_ODE ILUIAU 探索 : 一次元ではこの部分は手動 /*** *** FID_TS_ODE **/ statc vod FID_TS_ODE (nt pnt p2) { nt cou; for(=;<=ilu[p-];++){ f(p2 == IALU[p-][-]) return; cou=ilu[p-]+; IALU[p-][cou-]=p2; ILU[p-]=cou; return; 変数名サイズ内容 ILU IALU [] []{[LU] 各節点の非零非対角成分数 各節点の非零非対角成分 ( 列番号 )

92 FEM3D 92 節点探索 :FID_TS_ODE 一次元ではこの部分は手動 /*** *** FID_TS_ODE **/ statc vod FID_TS_ODE (nt pnt p2) { nt cou; for(=;<=ilu[p-];++){ f(p2 == IALU[p-][-]) return; 既に IALU に含まれている場合は 次のペアへ cou=ilu[p-]+; IALU[p-][cou-]=p2; ILU[p-]=cou; return;

93 FEM3D 93 節点探索 :FID_TS_ODE 一次元ではこの部分は手動 /*** *** FID_TS_ODE **/ statc vod FID_TS_ODE (nt pnt p2) { nt cou; for(=;<=ilu[p-];++){ f(p2 == IALU[p-][-]) return; cou=ilu[p-]+; IALU[p-][cou-]=p2; ILU[p-]=cou; return; IALU に含まれていない場合は ILU に を加えて IALU に格納

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

95 FEM3D 95 行列コネクティビティ生成 : MAT_CO0(/) FID_TS_ODE (nn); FID_TS_ODE (nn2); FID_TS_ODE (nn3); FID_TS_ODE (nn); FID_TS_ODE (nn5); FID_TS_ODE (nn6); FID_TS_ODE (nn7); for(n=0;n<;n++){ =ILU[n]; for (=0;<;++){ COL[]=IALU[n][]; msort(colcol2); for(=;>0;--){ IALU[n][-]= COL[COL2[-]-]; 各節点において IALU[][] が小さい番号から大きい番号に並ぶようにソート ( 単純なバブルソート ) せいぜい 00 程度のものをソートする

96 FEM3D 96 CRS 形式への変換 :MAT_CO #nclude <stdo.h> #nclude "pfem_utl.h" #nclude "allocate.h" etern FILE* fp_log; vod MAT_CO() { nt ; 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(=0;<ilu[];++){ =+ndelu[]; temlu[]=ialu[][]-; deallocate_vector(ilu); deallocate_vector(ialu); C nde[ ] nde[0] 0 FORTRA nde(0) 0 0 ILU[ ] nde( ) ILU( )

97 FEM3D 97 CRS 形式への変換 :MAT_CO #nclude <stdo.h> #nclude "pfem_utl.h" #nclude "allocate.h" etern FILE* fp_log; vod MAT_CO() { nt ; 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(=0;<ilu[];++){ =+ndelu[]; temlu[]=ialu[][]-; PLU=nde[] tem のサイズ非ゼロ非対角成分総数 deallocate_vector(ilu); deallocate_vector(ialu); MAT_CO

98 FEM3D 9 CRS 形式への変換 :MAT_CO #nclude <stdo.h> #nclude "pfem_utl.h" #nclude "allocate.h" etern FILE* fp_log; vod MAT_CO() { nt ; 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(=0;<ilu[];++){ =+ndelu[]; temlu[]=ialu[][]-; tem に 0 から始まる節点番号を記憶 deallocate_vector(ilu); deallocate_vector(ialu); MAT_CO

99 FEM3D 99 CRS 形式への変換 :MAT_CO #nclude <stdo.h> #nclude "pfem_utl.h" #nclude "allocate.h" etern FILE* fp_log; vod MAT_CO() { nt ; 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(=0;<ilu[];++){ =+ndelu[]; temlu[]=ialu[][]-; deallocate_vector(ilu); deallocate_vector(ialu); これらはもはや不要 MAT_CO

100 FEM3D 00 全体処理 /** 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() ;

101 FEM3D 0 MAT_ASS_MAI: 全体構成 do pn= 2 ガウス積分点番号 ( 方向 ) do jpn= 2 ガウス積分点番号 ( 方向 ) do pn= 2 ガウス積分点番号 ( 方向 ) ガウス積分点 ( 個 ) における形状関数 およびその 自然座標系 における微分の算出 enddo enddo enddo do cel= ICELTOT 要素ループ 節点の座標から ガウス積分点における 形状関数の 全体座標系 における微分 およびヤコビアンを算出 (JACOBI) do e= 局所節点番号 do je= 局所節点番号全体節点番号 :p jp A pjp のtemLUにおけるアドレス : j e do pn= 2 ガウス積分点番号 ( 方向 ) do jpn= 2 ガウス積分点番号 ( 方向 ) do pn= 2 ガウス積分点番号 ( 方向 ) 要素積分 要素行列成分計算 全体行列への足しこみ enddo enddo enddo enddo enddo enddo e

102 FEM3D 02 係数行列 :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 ; nt pjpp; nt pnjpnpn; nt cel; nt eje; nt SE; nt nn2n3nn5n6n7n; double SH; double QPQMEPEMTPTM; double XX2X3XX5X6X7X; double YY2Y3YY5Y6Y7Y; double ZZ2Z3ZZ5Z6Z7Z; 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]= e0; WEI[]= e0; POS[0]= e0; POS[]= e0;

103 FEM3D 03 係数行列 :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 ; nt pjpp; nt pnjpnpn; nt cel; nt eje; nt SE; nt nn2n3nn5n6n7n; double SH; double QPQMEPEMTPTM; double XX2X3XX5X6X7X; double YY2Y3YY5Y6Y7Y; double ZZ2Z3ZZ5Z6Z7Z; 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]= e0; WEI[]= e0; POS[0]= e0; POS[]= e0; POS: WEI: 積分点座標重み係数

104 FEM3D 0 係数行列 :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(p=0;p<2;p++){ QP=.e0 + POS[p]; QM=.e0 - POS[p]; EP=.e0 + POS[jp]; EM=.e0 - POS[jp]; TP=.e0 + POS[p]; TM=.e0 - POS[p]; SHAPE[p][jp][p][0]= Oth * QM * EM * TM; SHAPE[p][jp][p][]= Oth * QP * EM * TM; SHAPE[p][jp][p][2]= Oth * QP * EP * TM; SHAPE[p][jp][p][3]= Oth * QM * EP * TM; SHAPE[p][jp][p][]= Oth * QM * EM * TP; SHAPE[p][jp][p][5]= Oth * QP * EM * TP; SHAPE[p][jp][p][6]= Oth * QP * EP * TP; SHAPE[p][jp][p][7]= Oth * QM * EP * TP;

105 FEM3D 05 係数行列 :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(p=0;p<2;p++){ QP=.e0 + POS[p]; QM=.e0 - POS[p]; EP=.e0 + POS[jp]; EM=.e0 - POS[jp]; TP=.e0 + POS[p]; TM=.e0 - POS[p]; SHAPE[p][jp][p][0]= Oth * QM * EM * TM; SHAPE[p][jp][p][]= Oth * QP * EM * TM; SHAPE[p][jp][p][2]= Oth * QP * EP * TM; SHAPE[p][jp][p][3]= Oth * QM * EP * TM; SHAPE[p][jp][p][]= Oth * QM * EM * TP; SHAPE[p][jp][p][5]= Oth * QP * EM * TP; SHAPE[p][jp][p][6]= Oth * QP * EP * TP; SHAPE[p][jp][p][7]= Oth * QM * EP * TP; QP EP TP QM j EM j j TM

106 FEM3D 06 係数行列 :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(p=0;p<2;p++){ QP=.e0 + POS[p]; QM=.e0 - POS[p]; EP=.e0 + POS[jp]; EM=.e0 - POS[jp]; TP=.e0 + POS[p]; TM=.e0 - POS[p]; SHAPE[p][jp][p][0]= Oth * QM * EM * TM; SHAPE[p][jp][p][]= Oth * QP * EM * TM; SHAPE[p][jp][p][2]= Oth * QP * EP * TM; SHAPE[p][jp][p][3]= Oth * QM * EP * TM; SHAPE[p][jp][p][]= Oth * QM * EM * TP; SHAPE[p][jp][p][5]= Oth * QP * EM * TP; SHAPE[p][jp][p][6]= Oth * QP * EP * TP; SHAPE[p][jp][p][7]= Oth * QM * EP * TP;

107 FEM3D 07 係数行列 :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(p=0;p<2;p++){ QP=.e0 + POS[p]; QM=.e0 - POS[p]; EP=.e0 + POS[jp]; EM=.e0 - POS[jp]; TP=.e0 + POS[p]; TM=.e0 - POS[p]; SHAPE[p][jp][p][0]= Oth * QM * EM * TM; SHAPE[p][jp][p][]= Oth * QP * EM * TM; SHAPE[p][jp][p][2]= Oth * QP * EP * TM; SHAPE[p][jp][p][3]= Oth * QM * EP * TM; SHAPE[p][jp][p][]= Oth * QM * EM * TP; SHAPE[p][jp][p][5]= Oth * QP * EM * TP; SHAPE[p][jp][p][6]= Oth * QP * EP * TP; SHAPE[p][jp][p][7]= Oth * QM * EP * TP; ) ( ) ( ) ( ) ( 3 2 ) ( ) ( ) ( ) ( 7 6 5

108 FEM3D 0 係数行列 :MAT_ASS_MAI(3/6) PQ[jp][p][0]= - Oth * EM * TM; PQ[jp][p][]= + Oth * EM * TM; PQ[jp][p][2]= + Oth * EP * TM; PQ[jp][p][3]= - Oth * EP * TM; PQ[jp][p][]= - Oth * EM * TP; PQ[jp][p][5]= + Oth * EM * TP; PQ[jp][p][6]= + Oth * EP * TP; PQ[jp][p][7]= - Oth * EP * TP; PE[p][p][0]= - Oth * QM * TM; PE[p][p][]= - Oth * QP * TM; PE[p][p][2]= + Oth * QP * TM; PE[p][p][3]= + Oth * QM * TM; PE[p][p][]= - Oth * QM * TP; PE[p][p][5]= - Oth * QP * TP; PE[p][p][6]= + Oth * QP * TP; PE[p][p][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][]= + 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]; n=icelod[cel][3]; n5=icelod[cel][]; n6=icelod[cel][5]; n7=icelod[cel][6]; n=icelod[cel][7]; ( j ) l PQ( j ) ( j ) PE( ) l ( j ) PT ( l j) ( j ) ( j ) 2 ( j ) 3 ( j ) 3 ( j ) j j j j における形状関数の一階微分

109 FEM3D 09 係数行列 :MAT_ASS_MAI(3/6) PQ[jp][p][0]= - Oth * EM * TM; PQ[jp][p][]= + Oth * EM * TM; PQ[jp][p][2]= + Oth * EP * TM; PQ[jp][p][3]= - Oth * EP * TM; PQ[jp][p][]= - Oth * EM * TP; PQ[jp][p][5]= + Oth * EM * TP; PQ[jp][p][6]= + Oth * EP * TP; PQ[jp][p][7]= - Oth * EP * TP; PE[p][p][0]= - Oth * QM * TM; PE[p][p][]= - Oth * QP * TM; PE[p][p][2]= + Oth * QP * TM; PE[p][p][3]= + Oth * QM * TM; PE[p][p][]= - Oth * QM * TP; PE[p][p][5]= - Oth * QP * TP; PE[p][p][6]= + Oth * QP * TP; PE[p][p][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][]= + 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]; n=icelod[cel][3]; n5=icelod[cel][]; n6=icelod[cel][5]; n7=icelod[cel][6]; n=icelod[cel][7];

110 FEM3D 0 係数行列 :MAT_ASS_MAI(/6) nodlocal[0]= n; nodlocal[]= n2; nodlocal[2]= n3; nodlocal[3]= n; nodlocal[]= n5; nodlocal[5]= n6; nodlocal[6]= n7; nodlocal[7]= n; X=XYZ[n-][0]; X2=XYZ[n2-][0]; X3=XYZ[n3-][0]; X=XYZ[n-][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-][]; Y=XYZ[n-][]; Y5=XYZ[n5-][]; Y6=XYZ[n6-][]; Y7=XYZ[n7-][]; Y=XYZ[n-][]; 節点の節点番号 QVC= Oth*(X+X2+X3+X+X5+X6+X7+X+Y+Y2+Y3+Y+Y5+Y6+Y7+Y); Z=XYZ[n-][2]; Z2=XYZ[n2-][2]; Z3=XYZ[n3-][2]; Z=XYZ[n-][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 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[0]= n; nodlocal[]= n2; nodlocal[2]= n3; nodlocal[3]= n; nodlocal[]= n5; nodlocal[5]= n6; nodlocal[6]= n7; nodlocal[7]= n; X=XYZ[n-][0]; X2=XYZ[n2-][0]; X3=XYZ[n3-][0]; X=XYZ[n-][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-][]; Y=XYZ[n-][]; Y5=XYZ[n5-][]; Y6=XYZ[n6-][]; Y7=XYZ[n7-][]; Y=XYZ[n-][]; 節点の X 座標 節点の Y 座標 座標値 : 節点番号から 引く QVC= Oth*(X+X2+X3+X+X5+X6+X7+X+Y+Y2+Y3+Y+Y5+Y6+Y7+Y); Z=XYZ[n-][2]; Z2=XYZ[n2-][2]; Z3=XYZ[n3-][2]; Z=XYZ[n-][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 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[0]= n; nodlocal[]= n2; nodlocal[2]= n3; nodlocal[3]= n; nodlocal[]= n5; nodlocal[5]= n6; nodlocal[6]= n7; nodlocal[7]= n; X=XYZ[n-][0]; X2=XYZ[n2-][0]; X3=XYZ[n3-][0]; X=XYZ[n-][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-][]; Y=XYZ[n-][]; Y5=XYZ[n5-][]; Y6=XYZ[n6-][]; Y7=XYZ[n7-][]; Y=XYZ[n-][]; 節点の X 座標 節点の Y 座標 座標値 : 節点番号から 引く QVC= Oth*(X+X2+X3+X+X5+X6+X7+X+Y+Y2+Y3+Y+Y5+Y6+Y7+Y); Z=XYZ[n-][2]; Z2=XYZ[n2-][2]; Z3=XYZ[n3-][2]; Z=XYZ[n-][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 X X5 X6 X7 X Y Y2 Y3 Y Y5 Y6 Y7 Y Z Z2 Z3 Z Z5 Z6 Z7 Z); Q QVOL C C T Q 0 体積当たり発熱量は位置 ( メッシュの中心の座標 c c ) に依存

113 FEM3D 3 係数行列 :MAT_ASS_MAI(/6) nodlocal[0]= n; nodlocal[]= n2; nodlocal[2]= n3; nodlocal[3]= n; nodlocal[]= n5; nodlocal[5]= n6; nodlocal[6]= n7; nodlocal[7]= n; X=XYZ[n-][0]; X2=XYZ[n2-][0]; X3=XYZ[n3-][0]; X=XYZ[n-][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-][]; Y=XYZ[n-][]; Y5=XYZ[n5-][]; Y6=XYZ[n6-][]; Y7=XYZ[n7-][]; Y=XYZ[n-][]; 座標値 : 節点番号から 引く QVC= Oth*(X+X2+X3+X+X5+X6+X7+X+Y+Y2+Y3+Y+Y5+Y6+Y7+Y); Z=XYZ[n-][2]; Z2=XYZ[n2-][2]; Z3=XYZ[n3-][2]; Z=XYZ[n-][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 X X5 X6 X7 X Y Y2 Y3 Y Y5 Y6 Y7 Y Z Z2 Z3 Z Z5 Z6 Z7 Z); Q T QVOL C C T Q 0

114 FEM3D 係数行列 :MAT_ASS_MAI(/6) nodlocal[0]= n; nodlocal[]= n2; nodlocal[2]= n3; nodlocal[3]= n; nodlocal[]= n5; nodlocal[5]= n6; nodlocal[6]= n7; nodlocal[7]= n; X=XYZ[n-][0]; X2=XYZ[n2-][0]; X3=XYZ[n3-][0]; X=XYZ[n-][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-][]; Y=XYZ[n-][]; Y5=XYZ[n5-][]; Y6=XYZ[n6-][]; Y7=XYZ[n7-][]; Y=XYZ[n-][]; QVC= Oth*(X+X2+X3+X+X5+X6+X7+X+Y+Y2+Y3+Y+Y5+Y6+Y7+Y); Z=XYZ[n-][2]; Z2=XYZ[n2-][2]; Z3=XYZ[n3-][2]; Z=XYZ[n-][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 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(/) #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 XKREAL X5KREAL X6KREAL X7KREAL X KREAL YKREAL Y2KREAL Y3KREAL YKREAL Y5KREAL Y6KREAL Y7KREAL Y KREAL ZKREAL Z2KREAL Z3KREAL ZKREAL Z5KREAL Z6KREAL Z7KREAL Z) { /** calculates JACOBIA & IVERSE JACOBIA d/d d/d & d/d **/ nt pjpp; double dxdqdydqdzdqdxdedydedzdedxdtdydtdzdt; double coef; double aa2a3a2a22a23a3a32a33; for(p=0;p<2;p++){ for(jp=0;jp<2;jp++){ for(p=0;p<2;p++){ PX[p][jp][p][0]=0.0; PX[p][jp][p][]=0.0; PX[p][jp][p][2]=0.0; PX[p][jp][p][3]=0.0; PX[p][jp][p][]=0.0; PX[p][jp][p][5]=0.0; PX[p][jp][p][6]=0.0; PX[p][jp][p][7]=0.0; l l l 入力 l ~ 出力 l l l det J l l l 各ガウス積分点 [p][jp][p] における値

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

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

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

119 FEM3D 9 JACOBI(2/) PY[p][jp][p][0]=0.0; PY[p][jp][p][]=0.0; PY[p][jp][p][2]=0.0; PY[p][jp][p][3]=0.0; PY[p][jp][p][]=0.0; PY[p][jp][p][5]=0.0; PY[p][jp][p][6]=0.0; PY[p][jp][p][7]=0.0; PZ[p][jp][p][0]=0.0; PZ[p][jp][p][]=0.0; PZ[p][jp][p][2]=0.0; PZ[p][jp][p][3]=0.0; PZ[p][jp][p][]=0.0; PZ[p][jp][p][5]=0.0; PZ[p][jp][p][6]=0.0; PZ[p][jp][p][7]=0.0; J J J J 2 3 J J J J J J /** **/ DETERMIAT of the JACOBIA dxdq = PQ[jp][p][0]*X + PQ[jp][p][]*X2 + PQ[jp][p][2]*X3 + PQ[jp][p][3]*X + PQ[jp][p][]*X5 + PQ[jp][p][5]*X6 + PQ[jp][p][6]*X7 + PQ[jp][p][7]*X; dydq = PQ[jp][p][0]*Y + PQ[jp][p][]*Y2 + PQ[jp][p][2]*Y3 + PQ[jp][p][3]*Y + PQ[jp][p][]*Y5 + PQ[jp][p][5]*Y6 + PQ[jp][p][6]*Y7 + PQ[jp][p][7]*Y; dzdq = PQ[jp][p][0]*Z + PQ[jp][p][]*Z2 + PQ[jp][p][2]*Z3 + PQ[jp][p][3]*Z + PQ[jp][p][]*Z5 + PQ[jp][p][5]*Z6 + PQ[jp][p][6]*Z7 + PQ[jp][p][7]*Z; dxde = PE[p][p][0]*X + PE[p][p][]*X2 + PE[p][p][2]*X3 + PE[p][p][3]*X + PE[p][p][]*X5 + PE[p][p][5]*X6 + PE[p][p][6]*X7 + PE[p][p][7]*X; dxdq J dydq J dzdq J 2 3

120 FEM3D 20 /** IVERSE JACOBIA **/ JACOBI(3/) dyde = PE[p][p][0]*Y + PE[p][p][]*Y2 + PE[p][p][2]*Y3 + PE[p][p][3]*Y + PE[p][p][]*Y5 + PE[p][p][5]*Y6 + PE[p][p][6]*Y7 + PE[p][p][7]*Y; dzde = PE[p][p][0]*Z + PE[p][p][]*Z2 + PE[p][p][2]*Z3 + PE[p][p][3]*Z + PE[p][p][]*Z5 + PE[p][p][5]*Z6 + PE[p][p][6]*Z7 + PE[p][p][7]*Z; dxdt = PT[p][jp][0]*X + PT[p][jp][]*X2 + PT[p][jp][2]*X3 + PT[p][jp][3]*X + PT[p][jp][]*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]*Y + PT[p][jp][]*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]*Z + PT[p][jp][]*Z5 + PT[p][jp][5]*Z6 + PT[p][jp][6]*Z7 + PT[p][jp][7]*Z; DETJ[p][jp][p]= dxdq*(dyde*dzdt-dzde*dydt) + dydq*(dzde*dxdt-dxde*dzdt) + dzdq*(dxde*dydt-dyde*dxdt); coef=.0 / DETJ[p][jp][p]; 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][p]=fabs(DETJ[p][jp][p]); J J J J 2 3 J J J J J J

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

122 FEM3D 22 /** IVERSE JACOBIA JACOBI(3/) dyde = PE[p][p][0]*Y + PE[p][p][]*Y2 + PE[p][p][2]*Y3 + PE[p][p][3]*Y + PE[p][p][]*Y5 + PE[p][p][5]*Y6 + PE[p][p][6]*Y7 + PE[p][p][7]*Y; dzde = PE[p][p][0]*Z + PE[p][p][]*Z2 + PE[p][p][2]*Z3 + PE[p][p][3]*Z + PE[p][p][]*Z5 + PE[p][p][5]*Z6 + PE[p][p][6]*Z7 + PE[p][p][7]*Z; dxdt = PT[p][jp][0]*X + PT[p][jp][]*X2 + PT[p][jp][2]*X3 + PT[p][jp][3]*X + PT[p][jp][]*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]*Y + PT[p][jp][]*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]*Z + PT[p][jp][]*Z5 + PT[p][jp][5]*Z6 + PT[p][jp][6]*Z7 + PT[p][jp][7]*Z; DETJ[p][jp][p]= dxdq*(dyde*dzdt-dzde*dydt) + dydq*(dzde*dxdt-dxde*dzdt) + dzdq*(dxde*dydt-dyde*dxdt); **/ coef=.0 / DETJ[p][jp][p]; 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[p][jp][p]=fabs(DETJ[p][jp][p]);

123 FEM3D 23 JACOBI(/) /** set the d/dx d/dy & d/dz components **/ PX[p][jp][p][0]=a*PQ[jp][p][0]+a2*PE[p][p][0]+a3*PT[p][jp][0]; PX[p][jp][p][]=a*PQ[jp][p][]+a2*PE[p][p][]+a3*PT[p][jp][]; PX[p][jp][p][2]=a*PQ[jp][p][2]+a2*PE[p][p][2]+a3*PT[p][jp][2]; PX[p][jp][p][3]=a*PQ[jp][p][3]+a2*PE[p][p][3]+a3*PT[p][jp][3]; PX[p][jp][p][]=a*PQ[jp][p][]+a2*PE[p][p][]+a3*PT[p][jp][]; PX[p][jp][p][5]=a*PQ[jp][p][5]+a2*PE[p][p][5]+a3*PT[p][jp][5]; PX[p][jp][p][6]=a*PQ[jp][p][6]+a2*PE[p][p][6]+a3*PT[p][jp][6]; PX[p][jp][p][7]=a*PQ[jp][p][7]+a2*PE[p][p][7]+a3*PT[p][jp][7]; PY[p][jp][p][0]=a2*PQ[jp][p][0]+a22*PE[p][p][0]+a23*PT[p][jp][0]; PY[p][jp][p][]=a2*PQ[jp][p][]+a22*PE[p][p][]+a23*PT[p][jp][]; PY[p][jp][p][2]=a2*PQ[jp][p][2]+a22*PE[p][p][2]+a23*PT[p][jp][2]; PY[p][jp][p][3]=a2*PQ[jp][p][3]+a22*PE[p][p][3]+a23*PT[p][jp][3]; PY[p][jp][p][]=a2*PQ[jp][p][]+a22*PE[p][p][]+a23*PT[p][jp][]; PY[p][jp][p][5]=a2*PQ[jp][p][5]+a22*PE[p][p][5]+a23*PT[p][jp][5]; PY[p][jp][p][6]=a2*PQ[jp][p][6]+a22*PE[p][p][6]+a23*PT[p][jp][6]; PY[p][jp][p][7]=a2*PQ[jp][p][7]+a22*PE[p][p][7]+a23*PT[p][jp][7]; PZ[p][jp][p][0]=a3*PQ[jp][p][0]+a32*PE[p][p][0]+a33*PT[p][jp][0]; PZ[p][jp][p][]=a3*PQ[jp][p][]+a32*PE[p][p][]+a33*PT[p][jp][]; PZ[p][jp][p][2]=a3*PQ[jp][p][2]+a32*PE[p][p][2]+a33*PT[p][jp][2]; PZ[p][jp][p][3]=a3*PQ[jp][p][3]+a32*PE[p][p][3]+a33*PT[p][jp][3]; PZ[p][jp][p][]=a3*PQ[jp][p][]+a32*PE[p][p][]+a33*PT[p][jp][]; PZ[p][jp][p][5]=a3*PQ[jp][p][5]+a32*PE[p][p][5]+a33*PT[p][jp][5]; PZ[p][jp][p][6]=a3*PQ[jp][p][6]+a32*PE[p][p][6]+a33*PT[p][jp][6]; PZ[p][jp][p][7]=a3*PQ[jp][p][7]+a32*PE[p][p][7]+a33*PT[p][jp][7]; a a a a a a a a a

124 FEM3D 2 係数行列 :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]; =-; f( jp!= p ){ S=ndeLU[p-]; E=ndeLU[p ]; for( =S;<E;++){ f( temlu[] == jp- ){ =; brea; 全体行列の非対角成分 A p jp :tem におけるアドレス p= nodlocal[e] jp= nodlocal[je] から始まる節点番号 2

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

126 FEM3D 26 係数行列 :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]; =-; f( jp!= p ){ S=ndeLU[p-]; E=ndeLU[p ]; for( =S;<E;++){ f( temlu[] == jp- ){ =; brea; 要素マトリクス ( e ~j e ) 全体マトリクス ( p ~j p ) の関係 :temlu におけるアドレス

127 FEM3D 27 係数行列 :MAT_ASS_MAI(6/6) QV0= 0.e0; COEFj= 0.e0; for(pn=0;pn<2;pn++){ for(jpn=0;jpn<2;jpn++){ for(pn=0;pn<2;pn++){ coef= fabs(detj[pn][jpn][pn])*wei[pn]*wei[jpn]*wei[pn]; PX= PX[pn][jpn][pn][e]; PY= PY[pn][jpn][pn][e]; PZ= PZ[pn][jpn][pn][e]; PXj= PX[pn][jpn][pn][je]; PYj= PY[pn][jpn][pn][je]; PZj= PZ[pn][jpn][pn][je]; COEFj+= coef*cod0*(px*pxj+py*pyj+pz*pzj); SH= SHAPE[pn][jpn][pn][e]; QV0+= SH * QVOL * coef; f (jp==p) { D[p-]+= COEFj; B[p-]+= QV0*QVC; f (jp!= p) { AMAT[]+= COEFj; j j j det J ddd

128 FEM3D 2 係数行列 :MAT_ASS_MAI(6/6) QV0= 0.e0; COEFj= 0.e0; for(pn=0;pn<2;pn++){ for(jpn=0;jpn<2;jpn++){ for(pn=0;pn<2;pn++){ coef= fabs(detj[pn][jpn][pn])*wei[pn]*wei[jpn]*wei[pn]; PX= PX[pn][jpn][pn][e]; PY= PY[pn][jpn][pn][e]; PZ= PZ[pn][jpn][pn][e]; PXj= PX[pn][jpn][pn][je]; PYj= PY[pn][jpn][pn][je]; PZj= PZ[pn][jpn][pn][je]; COEFj+= coef*cod0*(px*pxj+py*pyj+pz*pzj); SH= SHAPE[pn][jpn][pn][e]; QV0+= SH * QVOL * coef; f (jp==p) { D[p-]+= COEFj; B[p-]+= QV0*QVC; f (jp!= p) { AMAT[]+= COEFj; j j M j L f W W W d d d f I ) ( ) ( d d d J j j j det

129 FEM3D 29 係数行列 :MAT_ASS_MAI(6/6) QV0= 0.e0; COEFj= 0.e0; for(pn=0;pn<2;pn++){ for(jpn=0;jpn<2;jpn++){ for(pn=0;pn<2;pn++){ coef= fabs(detj[pn][jpn][pn])*wei[pn]*wei[jpn]*wei[pn]; PX= PX[pn][jpn][pn][e]; PY= PY[pn][jpn][pn][e]; PZ= PZ[pn][jpn][pn][e]; PXj= PX[pn][jpn][pn][je]; PYj= PY[pn][jpn][pn][je]; PZj= PZ[pn][jpn][pn][je]; COEFj+= coef*cod0*(px*pxj+py*pyj+pz*pzj); SH= SHAPE[pn][jpn][pn][e]; QV0+= SH * QVOL * coef; f (jp==p) { D[p-]+= COEFj; B[p-]+= QV0*QVC; f (jp!= p) { AMAT[]+= COEFj; j j M j L f W W W d d d f I ) ( ) ( j j J W W W coef det d d d J j j j det

130 FEM3D 30 係数行列 :MAT_ASS_MAI(6/6) QV0= 0.e0; COEFj= 0.e0; for(pn=0;pn<2;pn++){ for(jpn=0;jpn<2;jpn++){ for(pn=0;pn<2;pn++){ coef= fabs(detj[pn][jpn][pn])*wei[pn]*wei[jpn]*wei[pn]; PX= PX[pn][jpn][pn][e]; PY= PY[pn][jpn][pn][e]; PZ= PZ[pn][jpn][pn][e]; PXj= PX[pn][jpn][pn][je]; PYj= PY[pn][jpn][pn][je]; PZj= PZ[pn][jpn][pn][je]; COEFj+= coef*cod0*(px*pxj+py*pyj+pz*pzj); SH= SHAPE[pn][jpn][pn][e]; QV0+= SH * QVOL * coef; j j j f (jp==p) { D[p-]+= COEFj; B[p-]+= QV0*QVC; f (jp!= p) { AMAT[]+= COEFj;

131 FEM3D 3 係数行列 :MAT_ASS_MAI(6/6) QV0= 0.e0; COEFj= 0.e0; for(pn=0;pn<2;pn++){ for(jpn=0;jpn<2;jpn++){ for(pn=0;pn<2;pn++){ coef= fabs(detj[pn][jpn][pn])*wei[pn]*wei[jpn]*wei[pn]; PX= PX[pn][jpn][pn][e]; PY= PY[pn][jpn][pn][e]; ( e) ( e) f ( e PZ= PZ[pn][jpn][pn][e]; PXj= PX[pn][jpn][pn][je]; PYj= PY[pn][jpn][pn][je]; PZj= PZ[pn][jpn][pn][je]; COEFj+= coef*cod0*(px*pxj+py*pyj+pz*pzj); SH= SHAPE[pn][jpn][pn][e]; QV0+= SH * QVOL * coef; f (jp==p) { D[p-]+= COEFj; B[p-]+= QV0*QVC; f (jp!= p) { AMAT[]+= COEFj; ) ( e) T f Q dv Q V QVOL C C QVC C C QV 0 QVOL f ( e) V QV 0QVC T dv

132 FEM3D 32 MAT_ASS_BC: 全体構成 do = 節点ループ ( ディリクレ ) 境界条件を設定する節点をマーク (IWKX) enddo do = 節点ループ f (IWKX().eq.) then マークされた節点だったら対応する右辺ベクトル (B) の成分 対角成分 (D) の成分の修正 ( 行 列 ) do = nde(-)+ nde() 対応する非零非対角成分 (AMAT) の成分の修正 ( 行 ) enddo endf enddo Z T=0@Z= ma do = 節点ループ do = ndelu (-)+ nde () f (IWKX(tem ()).eq.) then 対応する非零非対角成分の節点がマークされていたら対応する右辺ベクトル 非零非対角成分 (AMAT) の成分の修正 ( 列 ) endf enddo enddo X Y X Z Y

133 FEM3D 33 境界条件 :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 jnbb0cel; nt nn2n3nn5n6n7n; nt qq2q3qq5q6q7q; 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 ) brea; for( b=odgrp_idex[b0];b<odgrp_idex[b0+];b++){ n=odgrp_item[b]; IWKX[n-][0]=;

134 FEM3D 3 境界条件 :MAT_ASS_BC(2/2) for(n=0;n<;n++){ f( IWKX[n][0] == ){ B[n]= 0.e0; D[n]=.e0; for(=ndelu[n];<ndelu[n+];++){ AMAT[]= 0.e0; for(n=0;n<;n++){ for(=ndelu[n];<ndelu[n+];++){ f (IWKX[temLU[]][0] == ) { AMAT[]= 0.e0;

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

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

137 FEMD 37 プログラム :d.c(6/6) /* // // BOUDARY condtons // */ /* X=Xmn */ =0; js= Inde[]; AMat[jS]= 0.0; Dag[ ]=.0; Rhs [ ]= 0.0; T 0 =0 対角成分 = 右辺 =0 非対角成分 =0 for(=0;<plu;++){ f(item[]==0){amat[]=0.0; 復習 : 一次元問題

138 FEMD 3 プログラム :d.c(6/6) /* // // BOUDARY condtons // */ /* X=Xmn */ =0; js= Inde[]; AMat[jS]= 0.0; Dag[ ]=.0; Rhs [ ]= 0.0; T 0 =0 対角成分 = 右辺 =0 非対角成分 =0 ゼロクリア for(=0;<plu;++){ f(item[]==0){amat[]=0.0; 復習 : 一次元問題

139 FEMD 39 プログラム :d.c(6/6) /* // // BOUDARY condtons // */ /* X=Xmn */ =0; js= Inde[]; AMat[jS]= 0.0; Dag[ ]=.0; Rhs [ ]= 0.0; T 0 =0 対角成分 = 右辺 =0 非対角成分 =0 消去 ゼロクリア for(=0;<plu;++){ f(item[]==0){amat[]=0.0; 行列の対称性を保つため 第一種境界条件を適用している節点に対応する 列 を 右辺に移項して消去する ( 今の場合は非対角成分を 0 にするだけで良い ) 復習 : 一次元問題

140 FEMD 0 第一種境界条件が T 0 の場合 /* // // BOUDARY condtons // */ 行列の対称性を保つため 第一種境界条件を適用している節点に対応する 列 を 右辺に移項して消去する /* X=Xmn */ =0; js= Inde[]; AMat[jS]= 0.0; Dag[ ]=.0; Rhs [ ]= PHImn; Dag j j Inde[ j ] Inde[ j] Amat Item[ ] Rhs j for(j=;<;++){ for (=Inde[j];<Inde[j+];++){ f(item[]==0){ Rhs [j]= Rhs[j] Amat[]*PHImn; AMat[]= 0.0; 復習 : 一次元問題

141 FEMD 第一種境界条件が T 0 の場合 /* // // BOUDARY condtons // */ 行列の対称性を保つため 第一種境界条件を適用している節点に対応する 列 を 右辺に移項して消去する /* X=Xmn */ =0; js= Inde[]; AMat[jS]= 0.0; Dag[ ]=.0; Rhs [ ]= PHImn; Dag for(j=;<;++){ for (=Inde[j];<Inde[j+];++){ f(item[]==0){ Rhs [j]= Rhs[j] Amat[]*PHImn; AMat[]= 0.0; j Rhs Rhs j j j Inde[ j ] Inde[ j] Amat Amat s s Amat Item[ mn s s ] Item[ ] where Item[ s ] 0 復習 : 一次元問題

142 FEM3D 2 境界条件 :MAT_ASS_BC(2/2) for(n=0;n<;n++){ f( IWKX[n][0] == ){ B[n]= 0.e0; D[n]=.e0; for(=ndelu[n];<ndelu[n+];++){ AMAT[]= 0.e0; for(n=0;n<;n++){ for(=ndelu[n];<ndelu[n+];++){ f (IWKX[temLU[]][0] == ) { AMAT[]= 0.e0; IWKX[n-][0]= となる節点に対して対角成分 = 右辺 =0 非零対角成分 =0 ゼロクリア ここでやっていることも一次元の時と全く変わらない

143 FEM3D 3 境界条件 :MAT_ASS_BC(2/2) for(n=0;n<;n++){ f( IWKX[n][0] == ){ B[n]= 0.e0; D[n]=.e0; for(=ndelu[n];<ndelu[n+];++){ AMAT[]= 0.e0; IWKX[n-][0]= となる節点を非零非対角成分として有している節点に対して 右辺へ移項 当該非零非対角成分 =0 for(n=0;n<;n++){ for(=ndelu[n];<ndelu[n+];++){ f (IWKX[temLU[]][0] == ) { AMAT[]= 0.e0; 消去 ゼロクリア ここでやっていることも一次元の時と全く変わらない

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 **/ #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() ;

146 FEM3D 6 SOLVE #nclude <stdo.h> #nclude <strng.h> #nclude <math.h> #nclude "pfem_utl.h" #nclude "allocate.h" etern FILE *fp_log; etern vod CG(); vod SOLVE() { nt jl; nt ERROR ICFLAG=0; CHAR_LEGTH BUF; /** PARAMETERs **/ ITER = pfemiarra[0]; CG 法の最大反復回数 RESID = pfemrarra[0]; CG 法の収束判定値 /** ITERATIVE solver **/ CG (PLU D AMAT ndelu temlu B X RESID ITER &ERROR); ITERactual= ITER;

147 FEM3D 7 前処理付き共役勾配法 Precondtoned Conjugate Gradent Method (CG) Compute r (0) = b-[a] (0) for = 2 solve [M] (-) = r (-) - = r (-) (-) f = p () = (0) else 前処理 : 対角スケーリング end - = - / -2 p () = (-) + - endf q () = [A]p () = - /p () q () () = (-) + p () r () = r (-) - q () chec convergence r p (-)

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

149 FEM3D 9 #nclude <stdo.h> #nclude <math.h> #nclude "precson.h" #nclude "allocate.h" etern FILE *fp_log; CG ソルバー (/6) vod CG ( KIT KIT PLUKREAL D[] KREAL AMAT[]KIT ndelu[] KIT temlu[] KREAL B[]KREAL X[]KREAL RESIDKIT ITER KIT *ERROR) { nt j; nt elsleusu; double WVAL; double BRM20BRM2DRM20DRM2; double S_TIMEE_TIME; double ALPHABETA; double CC0RHORHO0RHO; nt terpre; KREAL **WW; KIT R=0Z=Q=P=2DD=3; KIT MAXIT; KREAL TOL;

150 FEM3D 50 #nclude <stdo.h> #nclude <math.h> #nclude "precson.h" #nclude "allocate.h" etern FILE *fp_log; CG ソルバー (/6) WW[][0]= WW[][R] {r WW[][]= WW[][Z] { WW[][]= WW[][Q] {q vod CG ( KIT KIT PLUKREAL D[] WW[][2]= WW[][P] {p KREAL AMAT[]KIT ndelu[] KIT temlu[] WW[][3]= KREAL B[]KREAL WW[][DD] X[]KREAL RESIDKIT /{D ITER KIT *ERROR) { nt j; nt elsleusu; double WVAL; double BRM20BRM2DRM20DRM2; double S_TIMEE_TIME; double ALPHABETA; double CC0RHORHO0RHO; nt terpre; KREAL **WW; KIT R=0Z=Q=P=2DD=3; KIT MAXIT; KREAL TOL; Compute r (0) = b-[a] (0) for = 2 solve [M] (-) = r (-) - = r (-) (-) f = end p () = (0) else - = - / -2 p () = (-) + - p (-) endf q () = [A]p () = - /p () q () () = (-) + p () r () = r (-) - q () chec convergence r