Microsoft PowerPoint - 3D-FEM-2.ppt [互換モード]

Transcription

1 三次元弾性解析コード (2/3) マトリクス生成 22 年夏学期 中島研吾 科学技術計算 Ⅰ(482-27) コンピュータ科学特別講義 Ⅰ(48-24)

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

3 FEM3D-Prt2 3 test メインプログラム nput_cntl 制御データ入力 nput_grd メッシュファイル入力 fnd_node 節点探索 mt_con 行列コネクティビティ生成 msort ソート mt_con 行列コネクティビティ生成 mt_ss_mnss mn 係数行列生成 cob ヤコビアン計算 mt_ss_bc 境界条件処理 三次元弾性解析コードfem3d の構成 solve33 線形ソルバー制御 recover_stress 応力計算 output_ucd 可視化処理 cg_3 CG 法計算 cob ヤコビアン計算

4 FEM3D-Prt2 4 fem3d: いくつかの特徴 非対角成分 上三角 下三角を分けて記憶 ndeltemlal ndeutemuau U ブロックとして記憶 ベクトル : 節点 3 成分 行列 : 各ブロック9 成分 行列の各成分ではなく 節点上の3 変数に基づくブロックとして処理する L

5 FEM3D-Prt2 5 ブロックとして記憶 (/3) 記憶容量が減る nde tem に関する記憶容量を削減できる

6 FEM3D-Prt2 6 ブロックとして記憶 (2/3) 計算効率 間接参照 ( メモリに負担 ) と計算の比が大きくなる ベクトル スカラー共に効く :2 倍以上の性能 連続領域 キャッシュに載る ループあたりの計算量増加 do 3* do Y() D()*X() X X(3*-2) do k nde(-) nde() X2 X(3*-) kk tem(k) X3 X(3*) enddo enddo Y() Y() AMAT(k)*X(kk) Y(3*-2) D(9*-8)*XD(9*-7)*X2D(9*-6)*X3 Y(3*-) D(9*-5)*XD(9*-4)*X2D(9*-3)*X3 Y(3*I ) D(9*-2)*XD(9*-)*X2D(9*I )*X3 do k nde(-) nde() enddo enddo kk tem(k) X X(3*kk-2) X2 X(3*kk-) X3 X(3*kk) Y(3*-2) Y(3*-2)AMAT(9*k-8)*XAMAT(9*k-7)*X2 & AMAT(9*k-6)*X3 Y(3*-) Y(3*-)AMAT(9*k-5)*XAMAT(9*k-4)*X2 & AMAT(9*k-3)*X3 Y(3*I ) Y(3*I )AMAT(9*k-2)*XAMAT(9*k-)*X2 & AMAT(9*k )*X3

7 FEM3D-Prt2 7 ブロックとして記憶 (3/3) 計算の安定化 対角成分で割るのではなく 対角ブロックの完全 LU 分解を求めて解く 特に悪条件問題で有効

8 FEM3D-Prt2 8 Globl 変数表 :pfem_utl.h/f(/3) 変数名種別サイズ I/O 内容 fnme C [8] I メッシュファイル名 P I I 節点数 ICELTOT I I 要素数 ODGRPtot I I 節点グループ数 XYZ R [][3] I 節点座標 ICELOD I [ICELTOT][8] I 要素コネクティビティ ODGRP_IDEX I [ODGRPtot] I 各節点グループに含まれる節点数 ( 累積 ) ODGRP_ITEM I [ODGRP_IDEX[ODG RPTOT]] I 節点グループに含まれる節点 ODGRP_AME C8 [ODGRP_IDEX[ODG RPTOT]] I 節点グループ名 L U I O 各節点非対角成分数 ( 上三角 下三角 ) PL PU I O 非対角成分総数 ( 上三角 下三角 ) D R [9*] O 全体行列 : 対角ブロック BX R [3*] O 右辺ベクトル 未知数ベクトル

9 FEM3D-Prt2 9 Globl 変数表 :pfem_utl.h/f(2/3) 変数名種別サイズ I/O 内容 ALUG R [9*] O 全体行列 : 対角ブロックの完全 LU 分解 ALAU R [9*PL][9*PU] O 全体行列 : 上 下三角ブロック成分 ndel ndeu I [] O 全体行列 : 非零非対角ブロック数 teml temu I [PL] [PU] O 全体行列 : 上 下三角ブロック ( 列番号 ) IL IU I [] O 各節点の上 下三角ブロック数 IALIAUIAU I [][L][][U][][U] O 各節点の上 下三角ブロック ( 列番号 ) IWKX I [][2] O ワーク用配列 METHOD I I 反復解法 ( に固定 ) PRECOD I I 前処理手法 (: ブロック SSOR: ブロック対角スケーリング ) ITER ITERctul I I 反復回数の上限 実際の反復回数 RESID R I 打ち切り誤差 (.e-8 に設定 ) SIGMA_DIAG R I LU 分解時の対角成分係数 (. に設定 ) pfemirr I [] O 諸定数 ( 整数 ) pfemrrr R [] O 諸定数 ( 実数 )

10 FEM3D-Prt2 Globl 変数表 :pfem_utl.h/f(3/3) 変数名種別サイズ I/O 内容 O8th R I.25 PQ PE PT R [2][2][8] O 各ガウス積分点における η ζ ( ~ 8) POS WEI R [2] O 各ガウス積分点の座標 重み係数 COL COL2 I [] O ソート用ワーク配列 SHAPE R [2][2][2][8] O 各ガウス積分点における形状関数 (~8) PX PY PZ R [2][2][2][8] O 各ガウス積分点における ( ~ 8 ) DETJ R [2][2][2] O 各ガウス積分点におけるヤコビアン行列式 ELAST POISSO R I ヤング率 ポアソン比 SIGMA_ TAU_ R [][3] O 節点における垂直 せん断応力成分

11 FEM3D-Prt2 用語の定義 ブロック ( 節点 ): 成分 ( 自由度 ):

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

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

14 FEM3D-Prt2 4 全体処理 #nclude <stdo.h> #nclude <stdlb.h> FILE* fp_log; #defne GLOBAL_VALUE_DEFIE #nclude "pfem_utl.h" etern vod IPUT_CTL(); etern vod IPUT_GRID(); etern vod MAT_CO(); etern vod MAT_CO(); etern vod MAT_ASS_MAI(); etern vod MAT_ASS_BC(); ASS etern vod SOLVE33(); etern vod RECOVER_STRESS(); etern vod OUTPUT_UCD(); nt mn() { /** Logfle for debug **/ f( (fp_logfopen("log.log""w")) ULL){ fprntf(stdout"nput fle cnnot be opened! n"); et(); IPUT_CTL(); IPUT_GRID(); MAT_CO(); MAT_CO(); MAT_ASS_MAI(); MAT_ASS_BC() ; SOLVE33(); RECOVER_STRESS(); OUTPUT_UCD() ; MAT_CO: ILIU IAL IAU 生成 MAT_CO: nde tem 生成 とりあえず から始まる節点番号を記憶

15 FEM3D-Prt2 5 MAT_CO: 全体構成 do cel ICELTOT 8 節点相互の関係から ILIUIALIAU を生成 (FID_ODE) enddo ( ) 8 7 ( ) ( ) 5 6 ) ( ) 4 3 ( ) 2 ( η ζ ) ( ) ( )

16 FEM3D-Prt2 6 行列コネクティビティ生成 : MAT_CO(/4) ( ) #nclude <stdo.h> #nclude "pfem_utl.h h" #nclude "llocte.h" etern FILE *fp_log; etern vod msort(nt* nt* nt); sttc vod FID_TS_ODE (ntnt); vod MAT_CO() { nt kceln; nt nn2n3n4n5n6n7n8; nt ; 2 256;( この変数は使用しない ) U 26; L 26; IL(KIT* )llocte_vector(seof(kit)); IAL(KIT**)llocte _ mtr(seof(kit)l); ( ) IU(KIT* )llocte_vector(seof(kit)); IAU(KIT**)llocte_mtr(seof(KIT)U); for(;<;) IL[]; for(;<;) for(;<l;) IAL[][]; for(;<;) IU[]; for(;<;) for(;<u;) IAU[][]; UL: 各節点における ( 上下 ) 非ゼロ非対角ブロックの最大数 ( 接続する節点数 ) 今の問題の場合はわかっているので このようにできる 不明の場合の実装 : レポート課題

17 FEM3D-Prt2 7 行列コネクティビティ生成 : MAT_CO(/4) ( ) #nclude <stdo.h> #nclude "pfem_utl.h h" #nclude "llocte.h" etern FILE *fp_log; etern vod msort(nt* nt* nt); sttc vod FID_TS_ODE (ntnt); vod MAT_CO() { nt kceln; nt nn2n3n4n5n6n7n8; nt ; 2 256;( この変数は使用しない ) U 26; L 26; 変数名サイズ内容各節点の上 下三角ブ IL IU [] ロック数 [][L] 各節点の上 下三角ブ IALIAU [][U] ロック ( 列番号 ) IL(KIT* )llocte_vector(seof(kit)); IAL(KIT**)llocte _ mtr(seof(kit)l); ( ) IU(KIT* )llocte_vector(seof(kit)); IAU(KIT**)llocte_mtr(seof(KIT)U); for(;<;) IL[]; for(;<;) for(;<l;) IAL[][]; for(;<;) IU[]; for(;<;) for(;<u;) IAU[][];

18 FEM3D-Prt2 8 行列コネクティビティ生成 : MAT_CO(2/4):( ) から始まる番号 for( cel;cel< ICELTOT;cel){ nicelod[cel][]; n2icelod[cel][]; n3icelod[cel][2]; n4icelod[cel][3]; n5icelod[cel][4]; n6icelod[cel][5]; n7icelod[cel][6]; n8icelod[cel][7]; FID_TS_ODE (nn2); FID_TS_ODE (nn3); FID_TS_ODE (nn4); FID_TS_ODE (nn5); FID_TS_ODE (nn6); FID_TS_ODE (nn7); FID_TS_ODE (nn8); ( ) 8 7 ( ) ( ) 5 6 ) ( ) 4 3 ( ) 2 ( η ζ ) ( ) ( ) FID_ TS_ ODE (n2n); FID_TS_ODE (n2n3); FID_TS_ODE (n2n4); FID_TS_ODE (n2n5); FID_TS_ODE (n2n6); FID_TS_ODE (n2n7); FID_TS_ODE ODE (n2n8); n8); FID_TS_ODE (n3n); FID_TS_ODE (n3n2); FID_TS_ODE (n3n4); FID_TS_ODE (n3n5); FID_TS_ODE ODE (n3n6); n6); FID_TS_ODE (n3n7); FID_TS_ODE (n3n8);

19 FEM3D-Prt2 9 節点探索 :FID_TS_ODE ILIUIALIAU 探索 : 一次元ではこの部分は手動 sttc vod FID_TS_ODE (nt pnt p2) { nt kkcou; f( p > p2 ){ for(kk;kk<il[p-];kk){ f(p2 IAL[p-][kk-]) return; couil[p-]; ] IAL[p-][cou-]p2; IL[p-]cou; return; f( p2 > p ) { for(kk;kk<iu[p-];kk){ f(p2 IAU[p-][kk-]) return; couiu[p-]; IAU[p-][cou-]p2; ][cou ] p2; IU[p-]cou; return; 変数名サイズ内容 IL IU [] 各節点の上 下三角ブロック数 IALIAU [][L] [][U] 各節点の上 下三角ブロック ( 列番号 )

20 FEM3D-Prt2 2 節点探索 :FID_TS_ODE 一次元ではこの部分は手動 sttc vod FID_TS_ODE (nt pnt p2) { nt kkcou; f( p > p2 ){ for(kk;kk<il[p-];kk){ f(p2 IAL[p-][kk-]) return; couil[p-]; ] IAL[p-][cou-]p2; IL[p-]cou; return; 既に IALIAUIAU に含まれている場合は 次のペアへ f( p2 > p ) { for(kk;kk<iu[p-];kk){ f(p2 IAU[p-][kk-]) return; couiu[p-]; IAU[p-][cou-]p2; ][cou ] p2; IU[p-]cou; 8 return;

21 FEM3D-Prt2 2 節点探索 :FID_TS_ODE 一次元ではこの部分は手動 sttc vod FID_TS_ODE (nt pnt p2) { nt kkcou; f( p > p2 ){ for(kk;kk<il[p-];kk){ f(p2 IAL[p-][kk-]) return; couil[p-]; ] IAL[p-][cou-]p2; IL[p-]cou; return; IALIAUIAU に含まれていない場合は IL IUにを加えて IALIAUに格納 f( p2 > p ) { for(kk;kk<iu[p-];kk){ f(p2 IAU[p-][kk-]) return; couiu[p-]; IAU[p-][cou-]p2; ][cou ] p2; IU[p-]cou; 8 return;

22 FEM3D-Prt2 22 行列コネクティビティ生成 : MAT_CO(3/4) ( ) FID_TS_ODE (n4n); FID_TS_ODE ODE (n4n2); FID_TS_ODE (n4n3); FID_TS_ODE (n4n5); FID_TS_ODE (n4n6); FID_TS_ODE (n4n7); FID_TS_ODE (n4n8); FID_TS_ODE (n5n); FID_TS_ODE (n5n2); FID_TS_ODE (n5n3); FID_TS_ODE (n5n4); FID_TS_ODE (n5n6); FID_TS_ODE ODE (n5n7); n7); FID_TS_ODE (n5n8); FID_TS_ODE (n6n); FID_TS_ODE (n6n2); FID_TS_ODE (n6n3); FID_TS_ODE ODE (n6n4); n4); FID_TS_ODE (n6n5); FID_TS_ODE (n6n7); FID_TS_ODE (n6n8); ( ) 8 7 ( ) ( ) 5 6 ) ( ) 4 3 ( ) 2 ( η ζ ) ( ) ( ) FID_TS_ODE (n7n); FID_TS_ODE ODE (n7n2); 72); FID_TS_ODE (n7n3); FID_TS_ODE (n7n4); FID_TS_ODE (n7n5); FID_TS_ODE (n7n6); FID_TS_ODE (n7n8);

23 FEM3D-Prt2 23 行列コネクティビティ生成 : MAT_CO(4/4) ( ) FID_TS_ODE (n8n); FID_TS_ODE ODE (n8n2); n2); FID_TS_ODE (n8n3); FID_TS_ODE (n8n4); FID_TS_ODE (n8n5); FID_TS_ODE (n8n6); FID_TS_ODE (n8n7); for(n;n<;n){ IL[n]; for (k;k<;k){ COL[k]IAL[n][k]; msort(colcol2); for(k;k>;k--){ IAL[n][-k] COL[COL2[k-]-]; 各節点において IAL[][k]IAU[][k] が小さい番号から大きい番号に並ぶようにソート ( 単純なバブルソート ) せいぜい 程度のものをソートする IU[n]; for (k;k<;k){ COL[k]IAU[n][k]; msort(colcol2); for(k;k>;k--){ IAU[n][-k] COL[COL2[k-]-];

24 FEM3D-Prt2 24 CRS 形式への変換 :MAT_CO #nclude <stdo.h> #nclude "pfem_utl.h" #nclude "llocte.h" etern FILE* fp_log; vod MAT_CO() { nt kkk; ndel(kit*)llocte_vector(seof(kit)); ndeu(kit*)llocte_vector(seof(kit)); ); for(;<;) ndel[]; for(;<;) ndeu[]; C ndel[ ] IL[ k] for(;<;){ ndel[]ndel[]il[]; k ndeu[]ndeu[]iu[]; PLndeL[]; ndeu[ ] PUndeU[]; IU[ k] k teml(kit*)llocte_vector(seof(kit)pl); temu(kit*)llocte_vector(seof(kit)pu); PU); ndel[] ndeu[] for(;<;){ for(k;k<il[];k){ kkkndel[]; teml[kk]ial[][k]; for(k;k<iu[];k){ kkkndeu[]; temu[kk]iau[][k]; dellocte_vector(il); dellocte_vector(iu); dellocte_vector(ial); dellocte_vector(iau); FORTRA ndel[ ] ndeu[ ] k k IL[ k ] IU[ k] ndel[] ndeu[]

25 FEM3D-Prt2 25 CRS 形式への変換 :MAT_CO #nclude <stdo.h> #nclude "pfem_utl.h" #nclude "llocte.h" etern FILE* fp_log; vod MAT_CO() { nt kkk; ndel(kit*)llocte_vector(seof(kit)); ndeu(kit*)llocte_vector(seof(kit)); ); for(;<;) ndel[]; for(;<;) ndeu[]; for(;<;){ ndel[]ndel[]il[]; ndeu[]ndeu[]iu[]; PLndeL[]; PUndeU[]; teml(kit*)llocte_vector(seof(kit)pl); temu(kit*)llocte_vector(seof(kit)pu); PU); for(;<;){ for(k;k<il[];k){ kkkndel[]; teml[kk]ial[][k]; for(k;k<iu[];k){ kkkndeu[]; temu[kk]iau[][k]; dellocte_vector(il); dellocte_vector(iu); dellocte_vector(ial); dellocte_vector(iau); PLndeL[] teml のサイズ非ゼロ非対角ブロック ( 下三角 ) 総数 PUndeU[] temuのサイズ非ゼロ非対角ブロック ( 上三角 ) 総数

26 FEM3D-Prt2 26 CRS 形式への変換 :MAT_CO #nclude <stdo.h> #nclude "pfem_utl.h" #nclude "llocte.h" etern FILE* fp_log; vod MAT_CO() { nt kkk; ndel(kit*)llocte_vector(seof(kit)); ndeu(kit*)llocte_vector(seof(kit)); ); for(;<;) ndel[]; for(;<;) ndeu[]; for(;<;){ ndel[]ndel[]il[]; ndeu[]ndeu[]iu[]; PLndeL[]; PUndeU[]; teml(kit*)llocte_vector(seof(kit)pl); temu(kit*)llocte_vector(seof(kit)pu); PU); for(;<;){ for(k;k<il[];k){ kkkndel[]; teml[kk]ial[][k]; for(k;k<iu[];k){ kkkndeu[]; temu[kk]iau[][k]; dellocte_vector(il); dellocte_vector(iu); dellocte_vector(ial); dellocte_vector(iau); temltemuにもから始まる節点番号を記憶

27 FEM3D-Prt2 27 CRS 形式への変換 :MAT_CO #nclude <stdo.h> #nclude "pfem_utl.h" #nclude "llocte.h" etern FILE* fp_log; vod MAT_CO() { nt kkk; ndel(kit*)llocte_vector(seof(kit)); ndeu(kit*)llocte_vector(seof(kit)); ); for(;<;) ndel[]; for(;<;) ndeu[]; for(;<;){ ndel[]ndel[]il[]; ndeu[]ndeu[]iu[]; PLndeL[]; PUndeU[]; teml(kit*)llocte_vector(seof(kit)pl); temu(kit*)llocte_vector(seof(kit)pu); PU); for(;<;){ for(k;k<il[];k){ kkkndel[]; teml[kk]ial[][k]; for(k;k<iu[];k){ kkkndeu[]; temu[kk]iau[][k]; dellocte_vector(il); dellocte_vector(iu); dellocte_vector(ial); dellocte_vector(iau); これらはもはや不要

28 FEM3D-Prt2 28 全体処理 #nclude <stdo.h> #nclude <stdlb.h> FILE* fp_log; #defne GLOBAL_VALUE_DEFIE #nclude "pfem_utl.h" etern vod IPUT_CTL(); etern vod IPUT_GRID(); etern vod MAT_CO(); etern vod MAT_CO(); etern vod MAT_ASS_MAI(); etern vod MAT_ASS_BC(); ASS etern vod SOLVE33(); etern vod RECOVER_STRESS(); etern vod OUTPUT_UCD(); nt mn() { /** Logfle for debug **/ f( (fp_logfopen("log.log""w")) ULL){ fprntf(stdout"nput fle cnnot be opened! n"); et(); IPUT_CTL(); IPUT_GRID(); MAT_CO(); MAT_CO(); MAT_ASS_MAI(); MAT_ASS_BC() ; SOLVE33(); RECOVER_STRESS(); OUTPUT_UCD() ;

29 FEM3D-Prt2 29 MAT_ ASS_ MAI: 全体構成 do kpn 2 ガウス積分点番号 (ζ 方向 ) do pn 2 ガウス積分点番号 (η 方向 ) do pn 2 ガウス積分点番号 ( 方向 ) ガウス積分点 (8 個 ) における形状関数 およびその 自然座標系 における微分の算出 enddo enddo enddo do cel ICELTOT 要素ループ 8 節点の座標から ガウス積分点における 形状関数の 全体座標系 における微分 およびヤコビアンを算出 (JACOBI) do e 8 局所節点番号 do e 8 局所節点番号全体節点番号 :p p A pp の temltemu におけるアドレス :kk do kpn 2 ガウス積分点番号 (ζ 方向 ) do pn 2 ガウス積分点番号 (η 方向 ) do pn 2 ガウス積分点番号 ( 方向 ) enddo enddo enddo enddo enddo enddo 要素積分 要素行列成分計算 全体行列への足しこみ e e

30 FEM3D-Prt2 3 ブロックとして記憶

31 FEM3D-Prt2 3 係数行列 :MAT_ASS_MAI(/7) #nclude <stdo.h> #nclude <mth.h> #nclude "pfem_utl.h" #nclude "llocte.h" etern FILE *fp_log; etern vod JACOBI(); vod MAT_ASS_MAI() { nt kkk; nt ppkp; nt pnpnkpn; nt cel; nt ee; nt SE; nt nn2n3n4n5n6n7n8; double vlavlbvlx; double QPQMEPEMTPTM; double XX2X3X4X5X6X7X8; double YY2Y3Y4Y5Y6Y7Y8; double ZZ2Z3Z4Z5Z6Z7Z8; double PXPYPZPXPYPZ; double VOL; double coef; double ; KIT nodlocal[8]; AL(KREAL*) llocte_vector(seof(kreal)9*pl); AU(KREAL*) llocte_vector(seof(kreal)9*pu); B (KREAL*) llocte_vector(seof(kreal)3* ); D (KREAL*) llocte_vector(seof(kreal)9*); X (KREAL*) llocte_vector(seof(kreal)3*); 下三角ブロック上三角ブロック右辺ベクトル対角ブロック未知数ベクトル for(;<9*pl;) AL[].; for(;<9*pu;) AU[].; ; for(;<3* ;) B[].; for(;<9* ;) D[].; for(;<3* ;) X[].;

32 FEM3D-Prt2 32 係数行列 :MAT_ASS_MAI(2/7) /*** ***/ WEI[].e; WEI[].e; POS[] e; POS[] e; POS: WEI: 積分点座標重み係数 IIT. PQ - st-order dervtve of shpe functon b QSI PE - st-order dervtve of shpe functon b ETA PT - st-order dervtve of shpe functon b ZET vla POISSO / (.e-poisso); vlb (.e-2.e*poisso)/(2.e*(.e-poisso)); vlx ELAST*(.e-POISSO)/((.ePOISSO)*(.e-2.e*POISSO)); vla vla * vlx; vlb vlb * vlx; for(p;p<2;p){ for(p;p<2;p){ for(kp;kp<2;kp){ QP.e POS[p]; QM.e - POS[p]; EP.e POS[p]; EM.e - POS[p]; TP.e POS[kp]; TM.e - POS[kp]; SHAPE[p][p][kp][] O8th * QM * EM * TM; SHAPE[p][p][kp][] O8th * QP * EM * TM; SHAPE[p][p][kp][2] O8th * QP * EP * TM; SHAPE[p][p][kp][3] O8th * QM * EP * TM; SHAPE[p][p][kp][4] O8th * QM * EM * TP; SHAPE[p][p][kp][5] O8th * QP * EM * TP; SHAPE[p][p][kp][6] O8th * QP * EP * TP; SHAPE[p][p][kp][7] O8th * QM * EP * TP;

33 FEM3D-Prt2 33 係数行列 :MAT_ASS_MAI(2/7) WEI[].e; WEI[].e; POS[] e; POS[] e; /*** ***/ IIT. PQ - st-order dervtve of shpe functon b QSI PE - st-order dervtve of shpe functon b ETA PT - st-order dervtve of shpe functon b ZET vla POISSO / (.e-poisso); vlb (.e-2.e*poisso)/(2.e*(.e-poisso)); vlx ELAST*(.e-POISSO)/((.ePOISSO)*(.e-2.e*POISSO)); vla vla * vlx; vlb vlb * vlx; for(p;p<2;p){ for(p;p<2;p){ for(kp;kp<2;kp){ QP.e POS[p]; QM.e - POS[p]; EP.e POS[p]; EM.e - POS[p]; TP.e POS[kp]; TM.e - POS[kp]; SHAPE[p][p][kp][] O8th * QM * EM * TM; SHAPE[p][p][kp][] O8th * QP * EM * TM; SHAPE[p][p][kp][2] O8th * QP * EP * TM; SHAPE[p][p][kp][3] O8th * QM * EP * TM; SHAPE[p][p][kp][4] O8th * QM * EM * TP; SHAPE[p][p][kp][5] O8th * QP * EM * TP; SHAPE[p][p][kp][6] O8th * QP * EP * TP; SHAPE[p][p][kp][7] O8th * QM * EP * TP;

34 FEM3D-Prt2 34 ひずみ 応力関係ず関係 ν ν ν ε ε ν ν ν ν ν ν σ σ ( )( ) ( ) ( ) E γ ε ν ν ν ν τ σ ( ) ( ) γ γ ν ν τ τ [ ] D { [ ]{ ε σ D

35 FEM3D-Prt2 35 係数行列 :MAT_ASS_MAI(2/7) WEI[].e; WEI[].e; POS[] e; POS[] e; /*** ***/ IIT. PQ - st-order dervtve of shpe functon b QSI PE - st-order dervtve of shpe functon b ETA PT - st-order dervtve of shpe functon b ZET vla POISSO / (.e-poisso); vlb (.e-2.e*poisso)/(2.e*(.e-poisso)); vlx ELAST*(.e-POISSO)/((.ePOISSO)*(.e-2.e*POISSO)); vla vla * vlx; vlb vlb * vlx; for(p;p<2;p){ for(p;p<2;p){ for(kp;kp<2;kp){ QP.e POS[p]; QM.e - POS[p]; EP.e POS[p]; EM.e - POS[p]; TP.e POS[kp]; TM.e - POS[kp]; vlx vla vlb SHAPE[p][p][kp][] O8th * QM * EM * TM; SHAPE[p][p][kp][] O8th * QP * EM * TM; SHAPE[p][p][kp][2] O8th * QP * EP * TM; SHAPE[p][p][kp][3] O8th * QM * EP * TM; SHAPE[p][p][kp][4] O8th * QM * EM * TP; SHAPE[p][p][kp][5] O8th * QP * EM * TP; SHAPE[p][p][kp][6] O8th * QP * EP * TP; SHAPE[p][p][kp][7] O8th * QM * EP * TP; E( ν ) ( ν )( 2ν ) ν E( ν ) ( ν ) ( ν )( 2ν ) ( 2ν ) E( ν ) 2( ν )( ν )( 2ν )

36 FEM3D-Prt2 36 係数行列 :MAT_ASS_MAI(2/7) WEI[].e; WEI[].e; POS[] e; POS[] e; /*** ***/ IIT. PQ - st-order dervtve of shpe functon b QSI PE - st-order dervtve of shpe functon b ETA PT - st-order dervtve of shpe functon b ZET vla POISSO / (.e-poisso); vlb (.e-2.e*poisso)/(2.e*(.e-poisso)); vlx ELAST*(.e-POISSO)/((.ePOISSO)*(.e-2.e*POISSO)); vla vla * vlx; vlb vlb * vlx; for(p;p<2;p){ for(p;p<2;p){ for(kp;kp<2;kp){ QP.e POS[p]; QM.e - POS[p]; EP.e POS[p]; EM.e - POS[p]; TP.e POS[kp]; TM.e - POS[kp]; vlx vla vlb SHAPE[p][p][kp][] O8th * QM * EM * TM; SHAPE[p][p][kp][] O8th * QP * EM * TM; SHAPE[p][p][kp][2] O8th * QP * EP * TM; SHAPE[p][p][kp][3] O8th * QM * EP * TM; SHAPE[p][p][kp][4] O8th * QM * EM * TP; SHAPE[p][p][kp][5] O8th * QP * EM * TP; SHAPE[p][p][kp][6] O8th * QP * EP * TP; SHAPE[p][p][kp][7] O8th * QM * EP * TP; E( ν ) ( ν )( 2ν ) ν E( ν ) Eν ( ν ) ( ν )( 2ν ) ( ν )( 2ν ) ( 2ν ) E( ν ) E 2( ν )( ν )( 2ν ) 2( ν )

37 FEM3D-Prt2 37 ひずみ 応力関係ず関係 ε ε σ σ vla vlx vla vla vla vlx γ ε τ σ lb vlb vlx vla vla γ γ τ τ vlb vlb [ ] D { [ ]{ ε σ D

38 FEM3D-Prt2 38 係数行列 :MAT_ASS_MAI(2/7) WEI[].e; WEI[].e; POS[] e; POS[] e; /*** ***/ IIT. PQ - st-order dervtve of shpe functon b QSI PE - st-order dervtve of shpe functon b ETA PT - st-order dervtve of shpe functon b ZET vla POISSO / (.e-poisso); vlb (.e-2.e*poisso)/(2.e*(.e-poisso)); vlx ELAST*(.e-POISSO)/((.ePOISSO)*(.e-2.e*POISSO)); vla vla * vlx; vlb vlb * vlx; for(p;p<2;p){ for(p;p<2;p){ for(kp;kp<2;kp){ QP.e POS[p]; QM.e - POS[p]; EP.e POS[p]; EM.e - POS[p]; TP.e POS[kp]; TM.e - POS[kp]; QP EP TP k () ( ) QM( ) ( ) () ( η ) EM() ( η ) ( ) ( ζ ) TM ( k ) ( ζ ) k ζ k SHAPE[p][p][kp][] O8th * QM * EM * TM; SHAPE[p][p][kp][] O8th * QP * EM * TM; SHAPE[p][p][kp][2] O8th * QP * EP * TM; SHAPE[p][p][kp][3] O8th * QM * EP * TM; SHAPE[p][p][kp][4] O8th * QM * EM * TP; SHAPE[p][p][kp][5] O8th * QP * EM * TP; SHAPE[p][p][kp][6] O8th * QP * EP * TP; SHAPE[p][p][kp][7] O8th * QM * EP * TP;

39 FEM3D-Prt2 39 係数行列 :MAT_ASS_MAI(2/7) WEI[].e; WEI[].e; POS[] e; POS[] e; /*** ***/ IIT. PQ - st-order dervtve of shpe functon b QSI PE - st-order dervtve of shpe functon b ETA PT - st-order dervtve of shpe functon b ZET vla POISSO / (.e-poisso); vlb (.e-2.e*poisso)/(2.e*(.e-poisso)); vlx ELAST*(.e-POISSO)/((.ePOISSO)*(.e-2.e*POISSO)); vla vla * vlx; vlb vlb * vlx; for(p;p<2;p){ for(p;p<2;p){ for(kp;kp<2;kp){ QP.e POS[p]; QM.e - POS[p]; EP.e POS[p]; EM.e - POS[p]; TP.e POS[kp]; TM.e - POS[kp]; SHAPE[p][p][kp][] O8th * QM * EM * TM; SHAPE[p][p][kp][] O8th * QP * EM * TM; SHAPE[p][p][kp][2] O8th * QP * EP * TM; SHAPE[p][p][kp][3] O8th * QM * EP * TM; SHAPE[p][p][kp][4] O8th * QM * EM * TP; SHAPE[p][p][kp][5] O8th * QP * EM * TP; SHAPE[p][p][kp][6] O8th * QP * EP * TP; SHAPE[p][p][kp][7] O8th * QM * EP * TP; ( ) 8 4 ( ) 7 3 ( ) ( ) 5 6 ( ) ( η ζ ) ( ) 2 ( ) ( )

40 FEM3D-Prt2 4 係数行列 :MAT_ASS_MAI(2/7) /*** ***/ WEI[].e; WEI[].e; POS[] e; POS[] e; IIT. PQ - st-order dervtve of shpe functon b QSI PE - st-order dervtve of shpe functon b ETA PT - st-order dervtve of shpe functon b ZET vla POISSO / (.e-poisso); 3 vlb (.e-2.e*poisso)/(2.e*(.e-poisso)); vlx ELAST*(.e-POISSO)/((.ePOISSO)*(.e-2.e*POISSO)); vla vla * vlx; vlb vlb * vlx; for(p;p<2;p){ for(p;p<2;p){ for(kp;kp<2;kp){ QP.e POS[p]; QM.e - POS[p]; EP.e POS[p]; EM.e - POS[p]; TP.e POS[kp]; TM.e - POS[kp]; SHAPE[p][p][kp][] O8th * QM * EM * TM; SHAPE[p][p][kp][] O8th * QP * EM * TM; SHAPE[p][p][kp][2] O8th * QP * EP * TM; SHAPE[p][p][kp][3] O8th * QM * EP * TM; SHAPE[p][p][kp][4] O8th * QM * EM * TP; SHAPE[p][p][kp][5] O8th * QP * EM * TP; SHAPE[p][p][kp][6] O8th * QP * EP * TP; SHAPE[p][p][kp][7] O8th * QM * EP * TP; 8 ( η ζ ) 8 ( η ζ ) 8 ( η ζ ) 8 ( η ζ ) 2 4 ( )( η)( ζ ) ( )( η)( ζ ) ( )( η )( ζ ) ( )( η)( ζ ) 5 ( η ζ ) 8 6( η ζ ) η ζ 8 7 ( η ζ ) η ζ 8 8( η ζ ) η ζ 8 ( )( η)( ζ ) ( )( )( ) ( )( )( ) ( )( )( )

41 FEM3D-Prt2 4 係数行列 :MAT_ASS_MAI(3/7) PQ[p][kp][] - O8th * EM * TM; PQ[p][kp][] O8th * EM * TM; PQ[p][kp][2] O8th * EP * TM; PQ[p][kp][3] - O8th * EP * TM; PQ[p][kp][4] - O8th * EM * TP; PQ[p][kp][5] O8th * EM * TP; PQ[p][kp][6] ][6] O8th * EP * TP; PQ[p][kp][7] - O8th * EP * TP; PE[p][kp][] - O8th * QM * TM; PE[p][kp][] - O8th * QP * TM; PE[p][kp][2] O8th * QP * TM; PE[p][kp][3] O8th * QM * TM; PE[p][kp][4] ][4] - O8th * QM * TP; PE[p][kp][5] - O8th * QP * TP; PE[p][kp][6] O8th * QP * TP; PE[p][kp][7] O8th * QM * TP; PT[p][p][] - O8th * QM * EM; [ p][p][ ] PT[p][p][2] - O8th * QP * EP; PT[p][p][3] - O8th * QM * EP; PT[p][p][4] O8th * QM * EM; PT[p][p][5] O8th * QP * EM; PT[p][p][6] O8th * QP * EP; PT[p][p][7] O8th * QM * EP; PQ( k) l ( η η ζ ζ k ) l PE ( k ) ( η η ζ ζ k ) η PT ( ) l ( η η ζ ζ k ) ζ ( ( )( ) PT[p][p][] - O8th * QP * EM; η ζ k ) η ζ k for( cel;cel< ICELTOT;cel){ nicelod[cel][]; n2icelod[cel][]; n3icelod[cel][2]; n4icelod[cel][3]; n5icelod[cel][4]; n6icelod[cel][5]; n7icelod[cel][6]; n8icelod[cel][7]; ( k η ζ ) 2 ( η ζ k ) ( η ζ k ) 8 3 ( η ζ k ) 8 ( η )( ζ ) ( η )( ζ ) ( η )( ζ ) における形状関数の一階微分 k k k

42 FEM3D-Prt2 42 係数行列 :MAT_ASS_MAI(3/7) PQ[p][kp][] - O8th * EM * TM; PQ[p][kp][] O8th * EM * TM; PQ[p][kp][2] O8th * EP * TM; PQ[p][kp][3] - O8th * EP * TM; PQ[p][kp][4] - O8th * EM * TP; PQ[p][kp][5] O8th * EM * TP; PQ[p][kp][6] ][6] O8th * EP * TP; PQ[p][kp][7] - O8th * EP * TP; PE[p][kp][] - O8th * QM * TM; PE[p][kp][] - O8th * QP * TM; PE[p][kp][2] O8th * QP * TM; PE[p][kp][3] O8th * QM * TM; PE[p][kp][4] ][4] - O8th * QM * TP; PE[p][kp][5] - O8th * QP * TP; PE[p][kp][6] O8th * QP * TP; PE[p][kp][7] O8th * QM * TP; PT[p][p][] - O8th * QM * EM; PT[p][p][] ] - O8th * QP * EM; PT[p][p][2] - O8th * QP * EP; PT[p][p][3] - O8th * QM * EP; PT[p][p][4] O8th * QM * EM; PT[p][p][5] O8th * QP * EM; PT[p][p][6] O8th * QP * EP; PT[p][p][7] O8th * QM * EP; for( cel;cel< ICELTOT;cel){ n ICELOD[cel][]; n2icelod[cel][]; n3icelod[cel][2]; n4icelod[cel][3]; n5icelod[cel][4]; n6icelod[cel][5]; n7icelod[cel][6]; n8icelod[cel][7]; ( ) ( ) nicelod[cel][]; ( ) 8 7 ( ) 5 6 ( ) 4 2 ( η ζ ) ( ) ( ) 3 ( )

43 FEM3D-Prt2 43 係数行列 :MAT_ASS_MAI(4/7) /** ** JACOBIA & IVERSE JACOBIA **/ nodlocal[] n; nodlocal[] n2; nodlocal[2] n3; nodlocal[3] n4; nodlocal[4] n5; nodlocal[5] n6; nodlocal[6] n7; nodlocal[7] n8; XXYZ[n-][]; XYZ[ ][] X2XYZ[n2-][]; X3XYZ[n3-][]; X4XYZ[n4-][]; X5XYZ[n5-][]; X6XYZ[n6-][]; ] X7XYZ[n7-][]; X8XYZ[n8-][]; 8 節点の節点番号 ( ) 8 7 ( ) ( ) 5 6 ( ) 8 節点の X 座標 ( ) ( η ζ ) ( ) ( ) ( ) YXYZ[n-][]; Y2XYZ[n2-][]; Y3XYZ[n3-][]; Y4XYZ[n4-][]; Y5XYZ[n5-][]; Y6XYZ[n6-][]; Y7XYZ[n7-][]; Y8XYZ[n8-][]; ZXYZ[n-][2]; Z2XYZ[n2-][2]; Z3XYZ[n3-][2]; Z4XYZ[n4-][2]; Z5XYZ[n5-][2]; Z6XYZ[n6-][2]; Z7XYZ[n7-][2]; Z8XYZ[n8-][2]; 8 節点の Y 座標 8 節点の Z 座標 座標値 : 節点番号から 引く JACOBI(DETJ PQ PE PT PX PY PZ X X2 X3 X4 X5 X6 X7 X8 Y Y2 Y3 Y4 Y5 Y6 Y7 Y8 Z Z2 Z3 Z4 Z5 Z6 Z7 Z8);

44 FEM3D-Prt2 44 #nclude <stdo.h> #nclude <mth.h> #nclude "precson.h" #nclude "llocte.h" /*** *** JACOBI ***/ JACOBI(/4) vod JACOBI( KREAL DETJ[2][2][2] KREAL PQ[2][2][8]KREAL PE[2][2][8]KREAL PT[2][2][8] KREAL PX[2][2][2][8]KREAL PY[2][2][2][8]KREAL PZ[2][2][2][8] KREAL XKREAL X2KREAL X3KREAL X4KREAL X5KREAL X6KREAL X7KREAL X8 { /** **/ KREAL YKREAL Y2KREAL Y3KREAL Y4KREAL Y5KREAL Y6KREAL Y7KREAL Y8 KREAL ZKREAL Z2KREAL Z3KREAL Z4KREAL Z5KREAL Z6KREAL Z7KREAL Z8) clcultes JACOBIA & IVERSE JACOBIA d/d d/d & d/d nt ppkp; double dxdqdydqdzdqdxdedydedzdedxdtdydtdzdt; double coef; double ; for(p;p<2;p){ for(p;p<2;p){ for(kp;kp<2;kp){ PX[p][p][kp][].; PX[p][p][kp][].; PX[p][p][kp][2].; ; PX[p][p][kp][3].; PX[p][p][kp][4].; PX[p][p][kp][5].; PX[p][p][kp][6].; PX[p][p][kp][7].; η ζ l l l 入力 ( )( l ~ 8) 出力 l l l det J l l l 各ガウス積分点 [p][p][kp] における値

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

46 FEM3D-Prt2 46 自然座標系における偏微分 (2/4) 自然座標系における偏微分 (2/4) マトリックス表示すると : [ ] J η η η η ζ ζ ζ ζ [ ] J : ヤコビのマトリクス J J J [ ] J : ヤコビのマトリクス (Jcob mtr Jcobn) [ ] J J J J J J J J J J ) J J J

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

48 FEM3D-Prt2 48 JACOBI(2/4) PY[p][p][kp][].; ; PY[p][p][kp][].; PY[p][p][kp][2].; PY[p][p][kp][3].; PY[p][p][kp][4].; PY[p][p][kp][5].; PY[p][p][kp][6].; ][k ][6] ; PY[p][p][kp][7].; ][k ][3] 2 3 PZ[p][p][kp][3].; [ J ] J J J PZ[p][p][kp][].; PZ[p][p][kp][].; PZ[p][p][kp][2].; PZ[p][p][kp][4].; PZ[p][p][kp][5].; PZ[p][p][kp][6].; PZ[p][p][kp][7].; J J 2 3 J J J J /** **/ DETERMIAT of the JACOBIA dxdq PQ[p][kp][]*X PQ[p][kp][]*X2 PQ[p][kp][2]*X3 PQ[p][kp][3]*X4 PQ[p][kp][4]*X5 PQ[p][kp][5]*X6 PQ[p][kp][6]*X7 PQ[p][kp][7]*X8; dydq PQ[p][kp][]*Y PQ[p][kp][]*Y2 PQ[p][kp][2]*Y3 PQ[p][kp][3]*Y4 PQ[p][kp][4]*Y5 PQ[p][kp][5]*Y6 PQ[p][kp][6]*Y7 PQ[p][kp][7]*Y8; dzdq PQ[p][kp][]*Z PQ[p][kp][]*Z2 PQ[p][kp][2]*Z3 PQ[p][kp][3]*Z4 PQ[p][kp][4]*Z5 PQ[p][kp][5]*Z6 PQ[p][kp][6]*Z7 PQ[p][kp][7]*Z8; dxde PE[p][kp][]*X PE[p][kp][]*X2 PE[p][kp][2]*X3 PE[p][kp][3]*X4 PE[p][kp][4]*X5 PE[p][kp][5]*X6 PE[p][kp][6]*X7 PE[p][kp][7]*X8; dxdq J dydq J 2 dzdq J 3

49 FEM3D-Prt2 49 JACOBI(3/4) dyde PE[p][kp][]*Y PE[p][kp][]*Y2 PE[p][kp][2]*Y3 PE[p][kp][3]*Y4 PE[p][kp][4]*Y5 PE[p][kp][5]*Y6 PE[p][kp][6]*Y7 PE[p][kp][7]*Y8; dzde PE[p][kp][]*Z PE[p][kp][]*Z2 PE[p][kp][2]*Z3 PE[p][kp][3]*Z4 PE[p][kp][4]*Z5 ][4] PE[p][kp][5]*Z6 ][5] PE[p][kp][6]*Z7 PE[p][kp][7]*Z8; dxdt PT[p][p][]*X PT[p][p][]*X2 PT[p][p][2]*X3 PT[p][p][3]*X4 PT[p][p][4]*X5 PT[p][p][5]*X6 PT[p][p][6]*X7 PT[p][p][7]*X8; dydt PT[p][p][]*Y ][] PT[p][p][]*Y2 ][] PT[p][p][2]*Y3 PT[p][p][3]*Y4 PT[p][p][4]*Y5 PT[p][p][5]*Y6 PT[p][p][6]*Y7 PT[p][p][7]*Y8; dzdt PT[p][p][]*Z PT[p][p][]*Z2 PT[p][p][2]*Z3 ] PT[p][p][3]*Z4 ] PT[p][p][4]*Z5 PT[p][p][5]*Z6 PT[p][p][6]*Z7 PT[p][p][7]*Z8; DETJ[p][p][kp] dxdq*(dyde*dzdt-dzde*dydt) dydq*(dzde*dxdt-dxde*dzdt) dzdq*(dxde*dydt-dyde*dxdt); J /** J2 J3 IVERSE JACOBIA [ J ] **/ J 2 J 22 J 23 coef. / DETJ[p][p][kp]; coef * ( dyde*dzdt - dzde*dydt ); J3 J32 J33 2 coef * ( dzdq*dydt - dydq*dzdt ); 3 coef * ( dydq*dzde - dzdq*dyde ); 2 coef * ( dzde*dxdt - dxde*dzdt ); 22 coef * ( dxdq*dzdt - dzdq*dxdt ); 23 coef * ( dzdq*dxde - dxdq*dzde ); 3 coef * ( dxde*dydt - dyde*dxdt ); 32 coef * ( dydq*dxdt - dxdq*dydt ); 33 coef * ( dxdq*dyde - dydq*dxde ); DETJ[p][p][kp]fbs(DETJ[p][p][kp]);

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

51 FEM3D-Prt2 5 /** **/ IVERSE JACOBIA JACOBI(3/4) dyde PE[p][kp][]*Y PE[p][kp][]*Y2 PE[p][kp][2]*Y3 PE[p][kp][3]*Y4 PE[p][kp][4]*Y5 PE[p][kp][5]*Y6 PE[p][kp][6]*Y7 PE[p][kp][7]*Y8; dzde PE[p][kp][]*Z PE[p][kp][]*Z2 PE[p][kp][2]*Z3 PE[p][kp][3]*Z4 PE[p][kp][4]*Z5 ][4] PE[p][kp][5]*Z6 ][5] PE[p][kp][6]*Z7 PE[p][kp][7]*Z8; dxdt PT[p][p][]*X PT[p][p][]*X2 PT[p][p][2]*X3 PT[p][p][3]*X4 PT[p][p][4]*X5 PT[p][p][5]*X6 PT[p][p][6]*X7 PT[p][p][7]*X8; dydt PT[p][p][]*Y ][] PT[p][p][]*Y2 ][] PT[p][p][2]*Y3 PT[p][p][3]*Y4 PT[p][p][4]*Y5 PT[p][p][5]*Y6 PT[p][p][6]*Y7 PT[p][p][7]*Y8; dzdt PT[p][p][]*Z PT[p][p][]*Z2 PT[p][p][2]*Z3 ] PT[p][p][3]*Z4 ] PT[p][p][4]*Z5 PT[p][p][5]*Z6 PT[p][p][6]*Z7 PT[p][p][7]*Z8; DETJ[p][p][kp] dxdq*(dyde*dzdt-dzde*dydt) dydq*(dzde*dxdt-dxde*dzdt) dzdq*(dxde*dydt-dyde*dxdt); coef. / DETJ[p][p][kp]; coef * ( dyde*dzdt - dzde*dydt ); 2-2 coef * ( dzdq*dydt dydq*dzdt ); 3 3 coef * ( dydq*dzde - dzdq*dyde ); [ J ] 2 coef * ( dzde*dxdt - dxde*dzdt ); coef * ( dxdq*dzdt - dzdq*dxdt ); 23 coef * ( dzdq*dxde - dxdq*dzde ); coef * ( dxde*dydt - dyde*dxdt ); 32 coef * ( dydq*dxdt - dxdq*dydt ); 33 coef * ( dxdq*dyde - dydq*dxde ); DETJ[p][p][kp]fbs(DETJ[p][p][kp]);

52 FEM3D-Prt2 52 /** set the d/dx d/dy & d/dz components JACOBI(4/4) **/ PX[p][p][kp][]*PQ[p][kp][]2*PE[p][kp][]3*PT[p][p][]; PX[p][p][kp][]*PQ[p][kp][]2*PE[p][kp][]3*PT[p][p][]; PX[p][p][kp][2]*PQ[p][kp][2]2*PE[p][kp][2]3*PT[p][p][2]; PX[p][p][kp][3]*PQ[p][kp][3]2*PE[p][kp][3]3*PT[p][p][3]; ][k ][3] PQ[ ][k ][3] 2 PE[ ][k ][3] 3 PT[ ][3]; PX[p][p][kp][4]*PQ[p][kp][4]2*PE[p][kp][4]3*PT[p][p][4]; PX[p][p][kp][5]*PQ[p][kp][5]2*PE[p][kp][5]3*PT[p][p][5]; PX[p][p][kp][6]*PQ[p][kp][6]2*PE[p][kp][6]3*PT[p][p][6]; PX[p][p][kp][7]*PQ[p][kp][7]2*PE[p][kp][7]3*PT[p][p][7]; PY[p][p][kp][]2*PQ[p][kp][]22*PE[p][kp][]23*PT[p][p][]; PY[p][p][kp][]2*PQ[p][kp][]22*PE[p][kp][]23*PT[p][p][]; ][k ][] PQ[ ][k ][] 22 PE[ ][k ][] 23 PT[ ][] PY[p][p][kp][2]2*PQ[p][kp][2]22*PE[p][kp][2]23*PT[p][p][2]; PY[p][p][kp][3]2*PQ[p][kp][3]22*PE[p][kp][3]23*PT[p][p][3]; PY[p][p][kp][4]2*PQ[p][kp][4]22*PE[p][kp][4]23*PT[p][p][4]; PY[p][p][kp][5]2*PQ[p][kp][5]22*PE[p][kp][5]23*PT[p][p][5]; PY[p][p][kp][6]2*PQ[p][kp][6]22*PE[p][kp][6]23*PT[p][p][6]; ] [p][ ] [ ] [ ] PY[p][p][kp][7]2*PQ[p][kp][7]22*PE[p][kp][7]23*PT[p][p][7]; PZ[p][p][kp][]3*PQ[p][kp][]32*PE[p][kp][]33*PT[p][p][]; PZ[p][p][kp][]3*PQ[p][kp][]32*PE[p][kp][]33*PT[p][p][]; PZ[p][p][kp][2]3*PQ[p][kp][2]32*PE[p][kp][2]33*PT[p][p][2]; PZ[p][p][kp][3]3*PQ[p][kp][3]32*PE[p][kp][3]33*PT[p][p][3]; PZ[p][p][kp][4]3*PQ[p][kp][4]32*PE[p][kp][4]33*PT[p][p][4]; PQ[p][kp][4] 32 PE[p][kp][4] 33 PT[p][p][4]; PZ[p][p][kp][5]3*PQ[p][kp][5]32*PE[p][kp][5]33*PT[p][p][5]; PZ[p][p][kp][6]3*PQ[p][kp][6]32*PE[p][kp][6]33*PT[p][p][6]; PZ[p][p][kp][7]3*PQ[p][kp][7]32*PE[p][kp][7]33*PT[p][p][7]; η ζ η ζ η ζ η ζ η ζ

53 FEM3D-Prt2 53 係数行列 :MAT_ASS_MAI(5/7) /** **/ COSTRUCT the GLOBAL MATRIX for(e;e<8;e){ pnodlocal[e]; for(e;e<8;e){ pnodlocal[e]; kk; f( p > p ){ SndeU[p-]; EndeU[p ]; for( ks;k<e;k){ f( temu[k-] p ){ kkk; brek; f( p < p ){ SndeL[p-]; EndeL[p ndel[p ]; for( ks;k<e;k){ f( teml[k-] p ){ kkk; brek; PX.e; PY.e; PZ.e; PX.e; PY.e; PZ.e; VOL.e; 全体行列の非対角ブロック A p p kk:temltemu におけるアドレス p nodlocal[e] p nodlocal[e] から始まる節点番号

54 FEM3D-Prt2 54 要素マトリクス :24 24 行列 各節点上の (uvw) 成分が物理的にも強くカップルしているので3 自由度をまとめて扱う :8 8 行列 e e ( ) 8 7 ( ) ( ) 5 6 ) ( ) 4 3 ( ) 2 ( η ζ ) ( ) ( ) e e e e e e 2 3 e e e e 2 22 e e 32 e e e e e e ( K8) e e

55 FEM3D-Prt2 55 全体マトリクス p p p p p p p p

56 FEM3D-Prt2 56 係数行列 :MAT_ASS_MAI(5/7) /** **/ COSTRUCT the GLOBAL MATRIX for(e;e<8;e){ pnodlocal[e]; for(e;e<8;e){ pnodlocal[e]; kk; f( p > p ){ SndeU[p-]; EndeU[p ]; for( ks;k<e;k){ f( temu[k-] p ){ kkk; brek; 要素マトリクス ( e ~ e ) 全体マトリクス ( p ~ p ) の関係 kk:temuteml におけるアドレスに を足した値 ( から始まる ) f( p < p ){ SndeL[p-]; EndeL[p ndel[p ]; for( ks;k<e;k){ f( teml[k-] p ){ kkk; brek; PX.e; PY.e; PZ.e; PX.e; PY.e; PZ.e; VOL.e;

57 FEM3D-Prt2 57 係数行列 :MAT_ASS_MAI(6/7).e; 2.e; 3.e; 2.e; 22.e; 23.e; 3.e; 32.e; 33.e; for(pn;pn<2;pn){ for(pn;pn<2;pn){ for(kpn;kpn<2;kpn){ coef -fbs(detj[pn][pn][kpn])*wei[pn]*wei[pn]*wei[kpn]; VOLcoef; PX PX[pn][pn][kpn][e]; PY PY[pn][pn][kpn][e]; PZ PZ[pn][pn][kpn][e]; PX PX[pn][pn][kpn][e]; PY PY[pn][pn][kpn][e]; PZ PZ[pn][pn][kpn][e]; (vlx*px*pxvlb*(py*pypz*pz))*coef; 22 (vlx*py*pyvlb*(pz*pzpx*px))*coef; 33 (vlx*pz*pzvlb*(px*pxpy*py))*coef; 2 (vla*px*py vlb*px*py)*coef; 3 (vla*px*pz vlb*px*pz)*coef; 2 (vla*py*px vlb*py*px)*coef; 23 (vla*py*pz vlb*py*pz)*coef; 3 (vla*pz*px vlb*pz*px)*coef; 32 (vla*pz*py PY vlb*pz*py)*coef; PY)

58 FEM3D-Prt2 58 係数行列 :MAT_ASS_MAI(6/7).e; 2.e; 3.e; 2.e; 22.e; 23.e; 3.e; 32.e; 33.e; for(pn;pn<2;pn){ ( ) J W W W ζ dt f for(pn ;pn<2;pn ){ for(pn;pn<2;pn){ for(kpn;kpn<2;kpn){ coef -fbs(detj[pn][pn][kpn])*wei[pn]*wei[pn]*wei[kpn]; VOLcoef; dv b D V ( ) k k J W W W ζ η coef det VOLcoef; PX PX[pn][pn][kpn][e]; PY PY[pn][pn][kpn][e]; PZ PZ[pn][pn][kpn][e]; PX PX[pn][pn][kpn][e]; PY PY[pn][pn][kpn][e]; ζ η d d d J b D ddd b D det PY PY[pn][pn][kpn][e]; PZ PZ[pn][pn][kpn][e]; (vlx*px*pxvlb*(py*pypz*pz))*coef; 22 (vlx*py*pyvlb*(pz*pzpx*px))*coef; 33 (vlx*pz*pzvlb*(px*pxpy*py))*coef; ζ η d d d J b D det 2 (vla*px*py vlb*px*py)*coef; 3 (vla*px*pz vlb*px*pz)*coef; 2 (vla*py*px vlb*py*px)*coef; 23 (vla*py*pz vlb*py*pz)*coef; 3 (vla*pz*px vlb*pz*px)*coef; 32 ( la PZ PY lb PZ PY) f ) ( ζ η f 32 (vla*pz*py vlb*pz*py)*coef; [ ] M L f W W W d d d f I ) ( ) ( ζ η ζ η ζ η [ ] k k k f W W W ) ( ζ η

59 FEM3D-Prt2 59 係数行列 :MAT_ASS_MAI(6/7) for(pn;pn<2;pn){ for(pn;pn<2;pn){ for(kpn;kpn<2;kpn){ coef -fbs(detj[pn][pn][kpn])*wei[pn]*wei[pn]*wei[kpn]; PX PX[pn][pn][kpn][e]; PY PY[pn][pn][kpn][e]; PZ PZ[pn][pn][kpn][e]; PX PX[pn][pn][kpn][e]; PY PY[pn][pn][kpn][e]; PZ PZ[pn][pn][kpn][e]; (vlx*px*pxvlb*(py*pypz*pz))*coef; 22 (vlx*py*pyvlb*(pz*pzpx*px))*coef; 33 (vlx*pz*pzvlb*(px*pxpy*py))*coef; σ vlx σ vla σ vla τ τ τ vla vlx vla 2 (vla*px*py vlb*px*py)*coef; 3 (vla*px*pz vlb*px*pz)*coef; 2 (vla*py*px vlb*py*px)*coef; 23 (vla*py*pz vlb*py*pz)*coef; 3 (vla*pz*px vlb*pz*px)*coef; 32 (vla*pz*py vlb*pz*py)*coef; vla vla vlx vlb vlb ε ε ε γ γ vlb γ

60 FEM3D-Prt2 6 X- 方向のつりあい式 (3/3) X 方向のつりあい式 (3/3) ( ) ν ν E b E ( )E ( )( ) ( )( ) ( ) ν ν ν ν ν ν ν E b E E D とすると [ ] [ ] [ ] [ ] [ ] { { { { { V b U b W V U D τ σ [ ] [ ] { { W b U b τ ( ) * ( ) [ ] [ ] [ ] [ ] [ ] [ ] ( ) { { * V T T T U dv b D [ ] [ ] [ ] [ ] ( ) { { [ ] [ ] [ ] [ ] { { V T T V T T V W dv b V dv b [ ] { V T dv X

61 FEM3D-Prt2 6 Y- 方向のつりあい式 [ ][ ] [ ] [ ] [ ] [ ] ( ) { { T T T T V T T T V dv b D [ ] [ ] [ ] [ ] ( ) { { [ ] [ ] [ ] [ ] { { [ ] { T V T T V T T dv Y W dv b U dv b [ ] { V dv Y Z- 方向のつりあい式 [ ][ ] [ ] [ ] [ ] [ ] ( ) { { T T T T V T T T V dv b U dv b W dv b D [ ] [ ] [ ] [ ] ( ) { { [ ] [ ] [ ] [ ] { { [ ] { T V T T V T T dv Z V dv b U dv b [ ] { V dv Z

62 FEM3D-Prt2 62 要素マトリクス :- 成分 X 方向 (/3) u v w vlxd; vla; vlbb ( K8) ( ) 8 7 ( ) 4 3 ( ) ( ) 5 6 ( ) 2 ( η ζ ) ( ) ( ) ( ) { vlx ( ) vlb V { vla vlb 2 { vla vlb 3 V V dv dv dv

63 FEM3D-Prt2 63 要素マトリクス :- 成分 X 方向 (/3) u v w vlxd; vla; vlbb ( K8) ( ) 8 7 ( ) 4 3 ( ) ( ) 5 6 ( ) 2 ( η ζ ) ( ) ( ) ( ) T [ ] [ ] [ ] T [ ] T b [ ] [ ] { D ( ) dv { U V T T [ ] [ ] b [ ] [ ] { ( ) dv { V V V T T [ ] [ ] b[ ] [ ] { dv{ W

64 FEM3D-Prt2 64 要素マトリクス :- 成分 Y 方向 (2/3) u v w vlxd; vla; vlbb ( K8) ( ) 8 7 ( ) 4 3 ( ) ( ) 5 6 ( ) 2 ( η ζ ) ( ) ( ) ( ) { vla vlb 2 22 V dv { vlx ( ) vlb { vla vlb 23 V V dv dv

65 FEM3D-Prt2 65 要素マトリクス :- 成分 Y 方向 (2/3) u v w vlxd; vla; vlbb ( K8) ( ) 8 7 ( ) 4 3 ( ) ( ) 5 6 ( ) 2 ( η ζ ) ( ) ( ) ( ) { [ ] [ ] ([ ] [ ]) T T b dv{ U V { D ( ) dv { V [ ] T [ T T ] b [ ] [ ] [ ] [ ] V b T T { [ ] [ ] [ ] [ ] dv { W V

66 FEM3D-Prt2 66 要素マトリクス :- 成分 Z 方向 (3/3) u v w vlxd; vla; vlbb ( K8) ( ) 8 7 ( ) 4 3 ( ) ( ) 5 6 ( ) 2 ( η ζ ) ( ) ( ) ( ) { vla vlb 3 V { vla vlb V dv dv { vld ( ) vlb V dv

67 FEM3D-Prt2 67 要素マトリクス :- 成分 Z 方向 (3/3) u v w vlxd; vla; vlbb ( K8) ( ) 8 7 ( ) 4 3 ( ) ( ) 5 6 ( ) 2 ( η ζ ) ( ) ( ) ( ) T T [ ] [ ] b [ ] [ ] V { ( ) dv{ U T T [ ] [ ] b[ ] [ ] { dv { V V V T T T [ ] [ ] b [ ] [ ] [ ] [ ] { D ( ) dv { W

68 FEM3D-Prt2 68 係数行列 :MAT_ASS_MAI(7/7) for(cel;cel<iceltot;cel){ for(e;e<8;e){ pnodlocal[e]; p p for(e;e<8;e){ <8 pnodlocal[e]; f (p > p) { AU[9*kk-9]; p p AU[9*kk-8]2; AU[9*kk-7]3; AU[9*kk-6]2; AU[9*kk-5]22; AU[9*kk-4]23; AU[9*kk-3]3; AU[9*kk-2]32; p p AU[9*kk-]33; f (p < p) { AL[9*kk-9]; AL[9*kk-8]2; AL[9*kk-7]3; AL[9*kk-6]2; p p AL[9*kk-5]22; AL[9*kk-4]23; AL[9*kk-3]3; AL[9*kk-2]32; AL[9*kk-]33; f (p p) { D[9*p-9]; 要素マトリクス D[9*p-8]2; ( e ~ e ) D[9*p-7]3; 全体マトリクス D[9*p-6]2; ( p ~ p ) の関係 D[9*p-5]22; D[9*p-4]23; D[9*p-3]3; kk:temuteml D[9*p-2]32; におけるアドレス に を足した値 D[9*p-]33; ( から始まる )

69 FEM3D-Prt2 69 要素マトリクス :24 24 行列 各節点上の (uvw) 成分が物理的にも強くカップルしているので3 自由度をまとめて扱う :8 8 行列 e e ( ) 8 7 ( ) ( ) 5 6 ) ( ) 4 3 ( ) 2 ( η ζ ) ( ) ( ) e e e e e e 2 3 e e e e 2 22 e e 32 e e e e e e ( K8) e e

70 FEM3D-Prt2 7 係数行列 :MAT_ASS_MAI(7/7) for(cel;cel<iceltot;cel){ for(e;e<8;e){ pnodlocal[e]; for(e;e<8;e){ <8 pnodlocal[e]; f (p > p) { AU[9*kk-9]; AU[9*kk-8]2; AU[9*kk-7]3; AU[9*kk-6]2; AU[9*kk-5]22; AU[9*kk-4]23; AU[9*kk-3]3; AU[9*kk-2]32; AU[9*kk-]33; f (p < p) { AL[9*kk-9]; AL[9*kk-8]2; AL[9*kk-7]3; AL[9*kk-6]2; AL[9*kk-5]22; AL[9*kk-4]23; AL[9*kk-3]3; AL[9*kk-2]32; AL[9*kk-]33; f (p p) { D[9*p-9]; D[9*p-8]2; D[9*p-7]3; D[9*p-6]2; D[9*p-5]22; D[9*p-4]23; D[9*p-3]3; D[9*p-2]32; D[9*p-]33;

71 FEM3D-Prt2 7 MAT_ ASS_ BC: 全体構成 do 節点ループ ( ディリクレ ) 境界条件を設定する節点をマーク (IWKX) enddo do 節点ループ f (IWKX().eq.) then マークされた節点だったら対応する右辺ベクトル (B) の成分 対角ブロック (D) の成分の修正 ( 行 列 ) do k ndel(-) ndel() 下三角ブロック対応する非対角 ( 下三角 ) ブロック (AL) の成分の修正 ( 行 ) enddo do k ndeu(-) ndeu() 上三角ブロック対応する非対角 ( 上三角 ) ブロック (AU) の成分の修正 ( 行 ) enddo endf enddo do 節点ループ do k ndel(-) ndel() 下三角ブロック f (IWKX(temL(k)).eq.) then 対応する非対角ブロックの節点がマークされていたら対応する右辺ベクトル 非対角 ( 下三角 ) ブロック (AL) の成分の修正 ( 列 ) endf enddo do k ndeu(-) ndeu() 上三角ブロック f (IWKX(temU(k)).eq.) then 対応する非対角ブロックの節点がマークされていたら対応する右辺ベクトル 非対角 ( 上三角 ) ブロック (AU) の成分の修正 ( 列 ) endf enddo enddo

72 FEM3D-Prt2 72 境界条件 :MAT_ASS_BC(/9) #nclude <stdo.h> #nclude <stdlb.h> #nclude <strng.h> #nclude "pfem_utl.h" #nclude "llocte.h" etern FILE *fp_log; vod MAT_ASS_BC() ASS { nt knbbcel; nt nn2n3n4n5n6n7n8; nt qq2q3q4q5q6q7q8; nt SE; double STRESSVAL; IWKX(KIT**) llocte_mtr(seof(kit)2); for(;<;) for(;<2;) IWKX[][];

73 FEM3D-Prt2 73 境界条件 :MAT_ASS_BC(2/9) /** **/ ZZm for(n;n<;n) IWKX[n][]; b-; for( b;b<odgrptot;b){ f( strcmp(odgrp_ame[b].nme"zm") ) brek; for( bodgrp_idex[b];b<odgrp_idex[b];b){ nodgrp_item[b]; IWKX[n-][]; for(n;n<;n){ f( IWKX[n][] ] ){ B[3*n ] B[3*n ] - D[9*n2]*.e; B[3*n] B[3*n] - D[9*n5]*.e; D[9*n2].e; D[9*n5].e; D[9*n6].e; D[9*n7].e; D[9*n8].e; B[3*n2].e; S ndel[n]; E ndel[n]; for(ks;k<e;k){ AL[9*k-3].e; AL[9*k-2].e; AL[9*k-].e; S ndeu[n]; E ndeu[n]; for(ks;k<e;k){ AU[9*k-3].e; AU[9*k-2].e; AU[9*k-].e; X Z Y U Z Z Y X

74 FEM3D-Prt2 74 境界条件 :MAT_ASS_BC(2/9) /** **/ ZZm for(n;n<;n) IWKX[n][]; b-; for( b;b<odgrptot;b){ f( strcmp(odgrp_ame[b].nme"zm") ) brek; for( bodgrp_idex[b];b<odgrp_idex[b];b){ nodgrp_item[b]; IWKX[n-][]; for(n;n<;n){ f( IWKX[n][] ] ){ B[3*n ] B[3*n ] - D[9*n2]*.e; B[3*n] B[3*n] - D[9*n5]*.e; D[9*n2].e; D[9*n5].e; D[9*n6].e; D[9*n7].e; D[9*n8].e; B[3*n2].e; S ndel[n]; E ndel[n]; for(ks;k<e;k){ S ndeu[n]; E ndeu[n]; AL[9*k-3].e; AL[9*k-2].e; AL[9*k-].e; for(ks;k<e;k){ AU[9*k-3].e; AU[9*k-2].e; AU[9*k-].e; 節点グループ名が Zm である節点において : IWKX[n-][] とする n:から始まる節点番号

75 FEM3D-Prt2 75 境界条件 :MAT_ASS_BC(2/9) /** **/ ZZm for(n;n<;n) IWKX[n][]; b-; for( b;b<odgrptot;b){ f( strcmp(odgrp_ame[b].nme"zm") ) brek; for( bodgrp_idex[b];b<odgrp_idex[b];b){ nodgrp_item[b]; IWKX[n-][]; for(n;n<;n){ f( IWKX[n][] ] ){ B[3*n ] B[3*n ] - D[9*n2]*.e; B[3*n] B[3*n] - D[9*n5]*.e; D[9*n2].e; D[9*n5].e; D[9*n6].e; D[9*n7].e; D[9*n8].e; B[3*n2].e; S ndel[n]; E ndel[n]; for(ks;k<e;k){ S ndeu[n]; E ndeu[n]; AL[9*k-3].e; AL[9*k-2].e; AL[9*k-].e; for(ks;k<e;k){ AU[9*k-3].e; AU[9*k-2].e; AU[9*k-].e; 節点グループ名が Zm である節点では Z 方向変位成分 にすなわち : B[3*n2]. (n ~-) FORTRA では B[3*n-2]. (n )

76 FEM3D-Prt2 76 境界条件 :MAT_ASS_BC(3/9) for(n;n<;n){ S ndel[n]; E ndel[n]; for(ks;k<e;k){ f (IWKX[temL[k-]-][] ) { B[3*n ] B[3*n ] - AL[9*k-7]*.e; B[3*n] B[3*n] - AL[9*k-4]*.e; B[3*n2] B[3*n2] - AL[9*k-]*.e; AL[9*k-7].e; AL[9*k-4].e; AL[9*k-].e; S ndeu[n]; E ndeu[n]; for(ks;k<e;k){ f (IWKX[temU[k-]-][] [ ][ ] ) { B[3*n ] B[3*n ] - AU[9*k-7]*.e; B[3*n] B[3*n] - AU[9*k-4]*.e; B[3*n2] B[3*n2] - AU[9*k-]*.e; AU[9*k-7].e; AU[9*k-4].e; AU[9*k-].e;

77 FEM3D-Prt2 77 一次元 で成立する方程式 u F ( mn ) m 方向にのみ自由度 ( 変位 u) 一様な : 断面積 A ヤング率 E 境界条件 :u( 固定 ) m : 大きさ F の力 ( 軸力 ) m 自重によるたわみ等はナシ : バネと同じ

78 FEM3D-Prt2 78 一次元 /* // プログラム :d.c(6/7) 境界条件 // BOUDARY condtons // */ u /* XXmn */ ; S Inde[]; AMt[S] ;.; Dg[ ].; Rhs [ ].; /* XXm */ / X Xm / -; Rhs[] F; for(k;k<plu;k){ f(item[k]){amt[k].; 対角成分 右辺 非対角成分

79 FEM3D-Prt2 79 一次元 /* // プログラム :d.c(6/7) 境界条件 // BOUDARY condtons // */ u /* XXmn */ ; S Inde[]; AMt[S] ;.; Dg[ ].; Rhs [ ].; /* XXm */ / X Xm / -; Rhs[] F; for(k;k<plu;k){ f(item[k]){amt[k].; 対角成分 右辺 非対角成分 ゼロクリア

80 FEM3D-Prt2 8 一次元 /* // プログラム :d.c(6/7) 境界条件 // BOUDARY condtons // */ u /* XXmn */ ; S Inde[]; AMt[S] ;.; Dg[ ].; Rhs [ ].; for(k;k<plu;k){ f(item[k]){amt[k].; ; /* XXm */ -; Rhs[] F; 行列の対称性を保つため 第一種境界条件を適用している節点に対応する 列 を 右辺に移項して消去する ( 今の場合は非対角成分をにするだけで良い ) 対角成分 右辺 非対角成分 消去 ゼロクリア

81 FEM3D-Prt2 8 一次元 /* // プログラム :d.c(6/7) 境界条件 // BOUDARY condtons // */ u /* XXmn */ ; S Inde[]; AMt[S] ;.; Dg[ ].; Rhs [ ].; for(k;k<plu;k){ f(item[k]){amt[k].; /* XXm */ -; Rhs[] F; 行列の対称性を保つため 第一種境界条件を適用している節点に対応する 列 を 右辺に移項して消去する ( 今の場合は非対角成分をにするだけで良い ) 対角成分 右辺 非対角成分 消去 ゼロクリア

82 FEM3D-Prt2 82 一次元 第一種境界条件が u の場合 /* // // BOUDARY condtons // */ /* XXmn */ ; S Inde[]; AMt[S].; Dg[ ].; Rhs [ ] Umn; 行列の対称性を保つため 第一種境界条件を適用している節点に対応する 列 を 右辺に移項して消去する for(;<;){ for (kinde[];k<inde[];k){ f(item[k]){ Rhs [] Rhs[] Amt[k]*Umn; AMt[k].; Inde[ ] Amt k u Item [ k ] k Inde[ ] Dg u Rhs

83 FEM3D-Prt2 83 一次元 第一種境界条件が u の場合 /* // // BOUDARY condtons // */ /* XXmn */ ; S Inde[]; AMt[S].; Dg[ ].; Rhs [ ] Umn; 行列の対称性を保つため 第一種境界条件を適用している節点に対応する 列 を 右辺に移項して消去する for(;<;){ for (kinde[];k<inde[];k){ f(item[k]){ Rhs [] Rhs[] Amt[k]*Umn; AMt[k].; Dg Rhs u Inde[ ] Amtk uitem [ k ] k Inde[ ] k k Amt k s u s Item[ k s ] Rhs Amt k s u mn where Item[ k s ]

84 FEM3D-Prt2 84 全体マトリクス p p p p p p p p

85 FEM3D-Prt2 85 同じことをやればよい 自分の 行 : 対角成分 それ以外 p p p p p p p p

86 FEM3D-Prt2 86 同じことをやればよい 自分の 列 : 右辺へ移項 非対角成分 p p p p p p p p

87 FEM3D-Prt2 87 境界条件 :MAT_ASS_BC(2/9) /** **/ ZZm for(n;n<;n) IWKX[n][]; b-; for( b;b<odgrptot;b){ f( strcmp(odgrp_ame[b].nme"zm") ) brek; for( bodgrp_idex[b];b<odgrp_idex[b];b){ nodgrp_item[b]; IWKX[n-][]; for(n;n<;n){ f( IWKX[n][] ] ){ B[3*n ] B[3*n ] - D[9*n2]*.e; B[3*n] B[3*n] - D[9*n5]*.e; D[9*n2].e; D[9*n5].e; D[9*n6].e; D[9*n7].e; D[9*n8].e; B[3*n2].e; S ndel[n]; E ndel[n]; for(ks;k<e;k){ S ndeu[n]; E ndeu[n]; AL[9*k-3].e; AL[9*k-2].e; AL[9*k-].e; for(ks;k<e;k){ AU[9*k-3].e; AU[9*k-2].e; AU[9*k-].e; 節点グループ名が Zm である節点において : IWKX[n-][] とする n:から始まる節点番号

88 FEM3D-Prt2 88 境界条件 :MAT_ASS_BC(2/9) /** **/ ZZm for(n;n<;n) IWKX[n][]; b-; for( b;b<odgrptot;b){ f( strcmp(odgrp_ame[b].nme"zm") ) brek; for( bodgrp_idex[b];b<odgrp_idex[b];b){ nodgrp_item[b]; IWKX[n-][]; for(n;n<;n){ f( IWKX[n][] ] ){ B[3*n ] B[3*n ] - D[9*n2]*.e; B[3*n] B[3*n] - D[9*n5]*.e; D[9*n2].e; D[9*n5].e; D[9*n6].e; D[9*n7].e; D[9*n8].e; B[3*n2].e; S ndel[n]; E ndel[n]; for(ks;k<e;k){ AL[9*k-3].e; AL[9*k-2].e; AL[9*k-].e; B[3*n] B[3*n] を変更してからD[9*n6] D[9*n7] をゼロクリア S ndeu[n]; E ndeu[n]; for(ks;k<e;k){ AU[9*k-3].e; AU[9*k-2].e; AU[9*k-].e; p p p p p p p p

89 FEM3D-Prt2 89 境界条件 :MAT_ASS_BC(2/9) /** **/ ZZm for(n;n<;n) IWKX[n][]; b-; for( b;b<odgrptot;b){ f( strcmp(odgrp_ame[b].nme"zm") ) brek; for( bodgrp_idex[b];b<odgrp_idex[b];b){ nodgrp_item[b]; IWKX[n-][]; p p p p for(n;n<;n){ f( IWKX[n][] ] ){ B[3*n ] B[3*n ] - D[9*n2]*.e; B[3*n] B[3*n] - D[9*n5]*.e; D[9*n2].e; D[9*n5].e; D[9*n6].e; D[9*n7].e; D[9*n8].e; B[3*n2].e; S ndel[n]; E ndel[n]; for(ks;k<e;k){ AL[9*k-3].e; AL[9*k-2].e; AL[9*k-].e; S ndeu[n]; E ndeu[n]; k strts from for(ks;k<e;k){ AU[9*k-3].e; AU[9*k-2].e; AU[9*k-].e; p p p p

90 FEM3D-Prt2 9 境界条件 :MAT_ASS_BC(3/9) for(n;n<;n){ S ndel[n]; E ndel[n]; for(ks;k<e;k){ f (IWKX[temL[k-]-][] ) { B[3*n ] B[3*n ] - AL[9*k-7]*.e; k strts from S ndeu[n]; E ndeu[n]; B[3*n] B[3*n] - AL[9*k-4]*.e; B[3*n2] B[3*n2] - AL[9*k-]*.e; AL[9*k-7].e; AL[9*k-4].e; AL[9*k-].e; 右辺を変更してからゼロクリア for(ks;k<e;k){ f (IWKX[temU[k-]-][] [ ][ ] ) { B[3*n ] B[3*n ] - AU[9*k-7]*.e; B[3*n] B[3*n] - AU[9*k-4]*.e; B[3*n2] B[3*n2] - AU[9*k-]*.e; AU[9*k-7].e; AU[9*k-4].e; AU[9*k-].e; p p p p p p

91 FEM3D-Prt2 9 w@zmn 自分の 行 : 対角成分 それ以外 p p p p p p p p

92 FEM3D-Prt2 92 境界条件 :MAT_ASS_BC(4/9) /** **/ ZZmn for(n;n<;n) IWKX[n][]; b-; for( b;b<odgrptot;b){ f( strcmp(odgrp_ame[b].nme"zmn") ) brek; w@zmn for( bodgrp_idex[b];b<odgrp_idex[b];b){ nodgrp_item[b]; IWKX[n-][]; for(n;n<;n){ f( IWKX[n][] ] ){ D[9*n2].e; D[9*n5].e; D[9*n6].e; D[9*n7].e; D[9*n8].e; B[3*n2].e; S ndel[n]; E ndel[n]; for(ks;k<e;k){ AL[9*k-3].e; AL[9*k-2].e; AL[9*k-].e; S ndeu[n]; E ndeu[n]; for(ks;k<e;k){ AU[9*k-3].e; AU[9*k-2].e; AU[9*k-].e; k strts from

93 FEM3D-Prt2 93 自分の 列 : 右辺へ移項 非対角成分 p p p p p p p p

94 FEM3D-Prt2 94 境界条件 :MAT_ASS_BC(5/9) for(n;n<;n){ S ndel[n]; E ndel[n]; for(ks;k<e;k){ f (IWKX[temL[k-]-][] ) { AL[9*k-7].e; AL[9*k-4].e; ; AL[9*k-].e; S ndeu[n]; E ndeu[n]; for(ks;k<e;k){ f (IWKX[temU[k-]-][] ) { AU[9*k-7].e; AU[9*k-4].e; AU[9*k-].e; w@zmn k strts t from

95 FEM3D-Prt2 95 自分の 行 : 対角成分 それ以外 p p p p p p p p

96 FEM3D-Prt2 96 境界条件 :MAT_ASS_BC(6/9) /** **/ XXmn for(n;n<;n) IWKX[n][]; b-; for( b;b<odgrptot;b){ f( strcmp(odgrp_ame[b].nme"xmn") ) brek; for( bodgrp_idex[b];b<odgrp_idex[b];b){ nodgrp_item[b]; IWKX[n-][]; for(n;n<;n){ f( IWKX[n][] ] ){ D[9*n ].e; D[9*n].e; D[9*n2].e; D[9*n3].e; D[9*n6].e; B[3*n ].e; u@xmn k strts from S ndel[n]; E ndel[n]; for(ks;k<e;k){ AL[9*k-9].e; AL[9*k-8].e; AL[9*k-7].e; S ndeu[n]; E ndeu[n]; for(ks;k<e;k){ AU[9*k-9].e; AU[9*k-8].e; AU[9*k-7].e;

97 FEM3D-Prt2 97 自分の 列 : 右辺へ移項 非対角成分 p p p p p p p p

98 FEM3D-Prt2 98 境界条件 :MAT_ASS_BC(7/9) for(n;n<;n){ S ndel[n]; E ndel[n]; for(ks;k<e;k){ f (IWKX[temL[k-]-][] ) { AL[9*k-9].e; AL[9*k-6].e; ; AL[9*k-3].e; u@xmn S ndeu[n]; E ndeu[n]; for(ks;k<e;k){ f (IWKX[temU[k-]-][] ) { AU[9*k-9].e; AU[9*k-6].e; AU[9*k-3].e; k strts from

99 FEM3D-Prt2 99 自分の 行 : 対角成分 それ以外 p p p p p p p p

100 FEM3D-Prt2 境界条件 :MAT_ASS_BC(8/9) /** **/ YYmn for(n;n<;n) IWKX[n][]; b-; for( b;b<odgrptot;b){ f( strcmp(odgrp_ame[b].nme"ymn") ) brek; v@ymn for( bodgrp_idex[b];b<odgrp_idex[b];b){ b<odgrp nodgrp_item[b]; IWKX[n-][]; for(n;n<;n){ f( IWKX[n][] ){ D[9*n].e; D[9*n4].e; D[9*n7].e; D[9*n3].e; D[9*n5].e; B[3*n].e; k strts from S ndel[n]; E ndel[n]; for(ks;k<e;k){ AL[9*k-6].e; AL[9*k-5].e; AL[9*k-4].e; S ndeu[n]; E ndeu[n]; for(ks;k<e;k){ AU[9*k-6].e; AU[9*k-5].e; AU[9*k-4].e;

101 FEM3D-Prt2 自分の 列 : 右辺へ移項 非対角成分 p p p p p p p p

102 FEM3D-Prt2 2 境界条件 :MAT_ASS_BC(9/9) for(n;n<;n){ S ndel[n]; E ndel[n]; for(ks;k<e;k){ f (IWKX[temL[k-]-][] ) { AL[9*k-8].e; AL[9*k-5].e; AL[9*k-2].e; S d ndeu[n]; U[] E ndeu[n]; for(ks;k<e;k){ f (IWKX[temU[k-]-][] ) { AU[9*k-8].e; AU[9*k-5].e; AU[9*k-2].e; v@ymn k strts from