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

FEM3D

2

対象とする問題：三次元定常熱伝導∂  ∂T  ∂  ∂T  ∂  ∂T  ɺ +  λλ +  λ + Q ( x, y , z ) = 0∂x  ∂x  ∂y  ∂y  ∂z  ∂z 

• 定常熱伝導＋発熱• 一様な熱伝導率 λ
• 直方体

Z

T=0@Z=zmax– 一辺長さ1の立方体（六面体）要素– 各方向にNX・NY・NZ個
NZ• 境界条件– T=0@Z=zmax

Y
NX
NY
X

• 体積当たり発熱量は位置（メッシュの中心の座標 xc,yc）に依存– Qɺ ( x, y, z ) = QVOL xC + yC

FEM3D

3

対象とする問題：三次元定常熱伝導∂  ∂T  ∂  ∂T  ∂  ∂T  ɺ +  λλ +  λ + Q ( x, y , z ) = 0∂x  ∂x  ∂y  ∂y  ∂z  ∂z 

• 原点から遠い部分が高温• 体積当たり発熱量は位置（メッシュの中心の座標）に依存– Qɺ ( x, y, z ) = xC + yC

movie

FEM3D

4

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

FEM3D

5

有限要素法の処理：プログラム• 初期化– 制御変数読み込み– 座標読み込み⇒要素生成（N:節点数，ICELTOT：要素数）– 配列初期化（全体マトリクス，要素マトリクス）– 要素⇒全体マトリクスマッピング（Index，Item）• マトリクス生成– 要素単位の処理（do icel= 1, ICELTOT）• 要素マトリクス計算• 全体マトリクスへの重ね合わせ– 境界条件の処理• 連立一次方程式– 共役勾配法（CG）

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

FEM3D7

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

FEM3D8

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

3

y
4

2
1
x

FEM3D9

自然/局所/要素 座標系（1/2）
Natural/Local/Element• 各要素を，自然座標系（ξ,η）の正方形要素[±1,±1] に変換する。
y

3

4

4

η

+1

-1

2

3

+1

ξ

1
x

1

-1

2

• 各要素の全体座標系（global coordinate）（x,y）における座標成分を，自然座標系における形状関数[N]（従属変数の内挿に使うのと同じ[N]）を使用して変換する

FEM3D10

自然/局所/要素 座標系（2/2）
Natural/Local/Element
y

3

4

4

η
+1

-1

2

3

+1

1
x

1

-1

• 各節点の座標：(x1,y1), (x2,y2), (x3,y3), (x4,y4)• 各節点における温度：T1, T2, T3, T4
4

T =  N i (ξ ,η ) ⋅ Ti
i =1
4

4

i =1

i =1

x =  N i (ξ ,η ) ⋅ xi , y =  N i (ξ ,η ) ⋅ yi2

ξ

FEM3D11

アイソパラメトリック要素
η

y

+1

-1
x

+1

ξ

-1

• 高次の補間関数を使えば，曲線，曲面も扱うことが可能となる。
• そういう意味で「自然座標系」と呼んでいる。
• Super-Parametric: Higher-Order Ni for (x,y)• Sub-Parametric: Lower-Order Ni for (x,y)

FEM3D

12

2D自然座標系の形状関数（1/3）η• 自然座標系における正方形上の内挿多項式は下式で与えら
れる：

4

+1

-1

T = α1 + α 2ξ + α 3η + α 4ξη• 各節点での条件より：
T1 + T2 + T3 + T4
− T1 + T2 + T3 − T4
α1 =, α2 =,
4
4− T1 − T2 + T3 + T4
T1 − T2 + T3 − T4α3 =, α4 =
4
4

1

3
+1

-1

2

ξ

FEM3D

13

2D自然座標系の形状関数（2/3）η

T = α1 + α 2ξ + α 3η + α 4ξη
T1 + T2 + T3 + T4 − T1 + T2 + T3 − T4
=
+
ξ+
4
4
− T1 − T2 + T3 + T4
T1 − T2 + T3 − T4
η+
ξη
4
4
1
1
= (1 − ξ − η + ξη )T1 + (1 + ξ − η − ξη )T2 +4
4
1(1 + ξ + η + ξη )T3 + 1 (1 − ξ + η − ξη )T4
4
4
1
1= (1 − ξ )(1 − η )T1 + (1 + ξ )(1 − η )T2 +4
4
1(1 + ξ )(1 + η )T3 + 1 (1 − ξ )(1 + η )T4
4
4

4

+1

-1
1

3
+1

-1

2

ξ

FEM3D

14

2D自然座標系の形状関数（2/3）η

T = α1 + α 2ξ + α 3η + α 4ξη
T1 + T2 + T3 + T4 − T1 + T2 + T3 − T4
=
+
ξ+
4
4
− T1 − T2 + T3 + T4
T1 − T2 + T3 − T4
η+
ξη
4
4
1
1
= (1 − ξ − η + ξη )T1 + (1 + ξ − η − ξη )T2 +4
4
1(1 + ξ + η + ξη )T3 + 1 (1 − ξ + η − ξη )T4
4
4
N1
N2
1
1
= (1 − ξ )(1 − η )T1 + (1 + ξ )(1 − η )T2 +4
4
1(1 + ξ )(1 + η )T3 + 1 (1 − ξ )(1 + η )T4
4
N3 4
N4

4

+1

-1
1

3
+1

-1

2

ξ

FEM3D

15

2D自然座標系の形状関数（3/3）η• 元の式に代入して，Tiについて整理すると以下のようになる：
T = N1T1 + N 2T2 + N 3T3 + N 4T4• 形状関数Niは以下のようになる：

4

+1

-1
1

1
(1 − ξ )(1 − η ), N 2 (ξ ,η ) = 1 (1 + ξ )(1 − η ),441
1
N 3 (ξ ,η ) = (1 + ξ )(1 + η ), N 4 (ξ ,η ) = (1 − ξ )(1 + η )4
4
N1 (ξ ,η ) =• 双一次（bi-linear）要素とも呼ばれる。
• 各節点におけるNiの値を計算してみよ

3
+1

-1

2

ξ

FEM3D16

三次元への拡張• 四面体要素：二次元における三角形要素– 任意の形状を扱うことができる。
– 特に一次要素は精度が悪く，一部の問題を除いてあまり使用されない。
– 高次の四面体要素は広く使用されている・・・• 本講義では低次六面体要素（アイソパラメトリック要素）を使用する。
– tri-linear• 自由度：温度@各節点上

FEM3D17

3D自然座標系の形状関数
1
N1 (ξ ,η , ζ ) = (1 − ξ )(1 − η )(1 − ζ )81
N 2 (ξ ,η , ζ ) = (1 + ξ )(1 − η )(1 − ζ )81
N 3 (ξ ,η , ζ ) = (1 + ξ )(1 + η )(1 − ζ )81
N 4 (ξ ,η , ζ ) = (1 − ξ )(1 + η )(1 − ζ )881
N 5 (ξ ,η , ζ ) = (1 − ξ )(1 − η )(1 + ζ )81
N 6 (ξ ,η , ζ ) = (1 + ξ )(1 − η )(1 + ζ )81
N 7 (ξ ,η , ζ ) = (1 + ξ )(1 + η )(1 + ζ )81
N 8 (ξ ,η , ζ ) = (1 − ξ )(1 + η )(1 + ζ )8(− 1,+1,+1) 8

T =  N i (ξ ,η , ζ ) ⋅ Ti
i =1(− 1,−1,+1) 58

x =  N i (ξ ,η , ζ ) ⋅ xi

7 (+ 1,+1,+1)
6 (+ 1,−1,+1)(− 1,+1,−1)

i =1
8

8

i =1

i =1

4

y =  N i (ξ ,η , ζ ) ⋅ yi , z =  N i (ξ ,η , ζ ) ⋅ zi
1

(ξ ,η , ζ ) = (− 1,−1,−1)

2

3 (+ 1,+1,−1)(+ 1,−1,−1)

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

FEM3D

19

ガラーキン法の適用（1/3）• 以下のような三次元熱伝導方程式を考慮する（熱伝導率一定）：
 ∂ 2T   ∂ 2T   ∂ 2T  ɺ λ 2  +  λ 2  +  λ 2  + Q = 0 ∂x   ∂y   ∂z 

T = [N ]{φ} （マトリクス形式），節点における温度を φ としてある。
要素内の温度分布• ガラーキン法に従い，重み関数を[N]とすると，各要素において以下の積分方程式が得られる：
  ∂ 2T   ∂ 2T   ∂ 2T  ɺ V [N ] λ  ∂x 2  + λ  ∂y 2  + λ  ∂z 2  + QdV = 0TFEM3D

20

ガラーキン法の適用（2/3）• 三次元のグリーンの定理
 ∂2B ∂2B ∂2B  ∂A ∂B ∂A ∂B ∂A ∂B 
∂B
 dV
+
+
A 2 + 2 + 2  dV = A
dS − ∂y∂z ∂n∂x ∂x ∂y ∂y ∂z ∂z 
 ∂x
V
S
V 







• 前式の2階微分の部分に適用（表面積分省略）：

 [N ] {λ (T ) + λ (T ) + λ (T )}dV
= −  {λ ([N ] T ) + λ ([N ] T ) + λ ([N ] T )}dV
T

, xx

, yy

, zz

V

T

,x

T

,x

,y

T

,y

,z

,z

V

• これに以下を代入する：[ ]T = [N ]{φ}, T, x = [N , x ]{φ}, T, y = N , y {φ}, T, z = [N , z ]{φ}FEM3D

21

ガラーキン法の適用（3/3）• 体積あたり発熱量の項 Qɺ を加えて次式が得られる：{() ([ ] [N ])+ λ ([N ] [N ])}dV ⋅ {φ}+  Qɺ [N ]dV = 0−  λ [N , x ] [N , x ] + λ N , y
V

T

T

T

,y

,z

,z

V

• この式を弱形式（weak form）と呼ぶ。元の微分方
程式では2階の微分が含まれていたが，上式では，グリーンの定理によって1階微分に低減されている。
– 弱形式によって近似関数（形状関数，内挿関数）に対する要求が弱くなっている：すなわち線形関数で2階微分の効果を記述できる。
– 項が増えただけで，一次元と同じ

FEM3D

22

境界条件を考慮した弱形式：各要素[k ] {φ }
(e)

[k ]

(e)

(e)

= {f }(e)()+  λ ([N ] [N ])dV
V

T

,z

V

[f ]

(e)

([ ] [N ])dV

=  λ [N , x ] [N , x ] dV +  λ N , yTTɺ=  Q[N ] dV
V

,z

V

T

,y

FEM3D

23

要素マトリクス：8×8行列
j[k ] (i, j = 1…8)ij

i(− 1,+1,+1) 8(− 1,−1,+1) 5

7 (+ 1,+1,+1)
6 (+ 1,−1,+1)(− 1,+1,−1)
4
1(ξ ,η , ζ ) = (− 1,−1,−1)

2

3 (+ 1,+1,−1)(+ 1,−1,−1)

FEM3D24

要素マトリクス：kij[k ] (i, j = 1…8)(− 1,+1,+1) 8

ij

7 (+ 1,+1,+1)(− 1,−1,+1) 5

6 (+ 1,−1,+1)(− 1,+1,−1)
4
1(ξ ,η , ζ ) = (− 1,−1,−1)[k ](e) =  λ ([N , x ]T [N , x ])dV +  λ ([N , y ]T [N , y ])dVV()+  λ [N , z ] [N , z ] dVTVVkij = −  {λ ⋅ N i , x ⋅ N j , x + λ ⋅ N i , y ⋅ N j , y + λ ⋅ N i , z ⋅ N j , z }dVV2

3 (+ 1,+1,−1)(+ 1,−1,−1)1D-Part225

数値的に積分を実施する方法• 台形公式• シンプソンの公式• ガウスの積分公式（ガウス＝ルジャンドル（Gauss-Legendre）とも呼ばれる，精度良い）• 不定積分を解析的に求めるのではなく，有限な数のサンプル点における値を利用する
X2m

 f (x )dx   [w ⋅ f (x )]

X1

k =1

k

k

1D-Part2

26

ガウス積分公式：一次元の例
シンプソンの公式より精度良い
Simpson’sGaussf ( x) = sin( x)X1

X 1 = 0,

X2

X2 =

π
4

-1 ξ 0 ξ +1

X3

,

1

X3 =

h = X 2 − X1 = X 3 − X1 =

π
2

π
4

h
S = [ f ( X 1 ) + 4 f ( X 2 ) + f ( X 3 )] = 1.0023
32

ξ1 , ξ 2 = ±0.5773502692
π /2

+1

0

−1

S =  f ( x )dx =  f (ξ ) h dξ2≅ h Wk ⋅ f (ξ k ) = 0.99847
k =1

1D-Part2

27

ガウスの積分公式• 無次元化された自然座標系 [-1,+1] を対象とする。
• m個の積分点を使用すると（2m-1）次の関数までは近似可能（従って1次，2次の内挿関数（形状関数）を使用するときは，m=2で十分）+1m f (ξ )dξ =  [w ⋅ f (ξ )]

−1

k =1

k

ξ=-1

ξ=0

ξ=+1

k

m = 1 ξ k = 0.00, wk = 2.00
m = 2 ξ k = ±0.577350, wk = 1.00
m = 3 ξ k = 0.00, wk = 8 / 9

ξ k = ±0.774597, wk = 5 / 9

ξｋ=-0.577350 ξｋ=+0.577350

FEM3D

28

Gaussian Quadrature
ガウスの積分公式• 自然座標系（ξ,η,ζ）上で定義⇒ガウス積分公式が使える（三次元）+1 +1 +1

I=

  f (ξ ,η , ζ ) dξdηdζ

−1 −1 −1

L

M

N

i =1

j =1

k =1

[

=    Wi ⋅W j ⋅ Wk ⋅ f (ξ i ,η j , ζ k )L,M,N： ξ,η,ζ 方向の積分点数(ξ i ,η j , ζ k ) ： 積分点の座標値
Wi , W j , Wk ： 積分点での重み係数]

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

30

Gaussian Quadrature
ガウスの積分公式
この組み合わせがよく使われる，二次元では4点におけるfの値を計算して足せば良い（三次元では8点）+1 +1I=  f (ξ ,η ) dξdη =   [W ⋅W ⋅ f (ξ ,η )]
m

n

i

−1 −1

i =1

j

i

j

j =1

= 1.0 ×1.0 × f (−0.57735,−0.57735) + 1.0 ×1.0 × f (−0.57735,+0.57735)+ 1.0 × 1.0 × f (+0.57735,+0.57735) + 1.0 × 1.0 × f (+0.57735,−0.57735)

FEM3D

31

あとは積分を求めれば良い• 自然座標系（ξ,η,ζ）上で定義⇒ガウス積分公式が使える（三次元）・・・しかし，微分が+1 +1 +1

I=

  f (ξ ,η , ζ ) dξdηdζ

−1 −1 −1

L

M

N

i =1

j =1

k =1

[

=    Wi ⋅W j ⋅ Wk ⋅ f (ξ i ,η j , ζ k )L,M,N： ξ,η,ζ 方向の積分点数(ξ i ,η j , ζ k ) ： 積分点の座標値
Wi , W j , Wk ： 積分点での重み係数]

FEM3D32

自然座標系における偏微分（1/4）• 偏微分の公式より以下のようになる：∂N i (ξ ,η , ζ ) ∂N i ∂x ∂N i ∂y ∂N i ∂z=++
∂x ∂ξ ∂y ∂ξ ∂z ∂ξ∂ξ∂N i (ξ ,η , ζ ) ∂N i ∂x ∂N i ∂y ∂N i ∂z=++∂η∂x ∂η ∂y ∂η ∂z ∂η
∂N i (ξ ,η , ζ ) ∂N i ∂x ∂N i ∂y ∂N i ∂z
=
+
+
∂ζ
∂x ∂ζ
∂y ∂ζ
∂z ∂ζ
 ∂N i ∂N i ∂N i  ∂ξ , ∂η , ∂ζ  は定義より簡単に求められるが
 ∂N i ∂N i ∂N i  ∂x , ∂y , ∂z  を実際の計算で使用する


FEM3D

33

自然座標系における偏微分（2/4）• マトリックス表示すると： ∂N i   ∂x ∂ξ   ∂ξ  ∂N i   ∂x= ∂η   ∂η ∂N i   ∂x ∂ζ   ∂ζ

∂y
∂ξ
∂y
∂η
∂y
∂ζ

∂z   ∂N i  ∂N i  ∂x ∂ξ   ∂x  ∂N ∂z   ∂N i  i
= [J ]∂η  ∂y ∂y 
∂z   ∂N i ∂N i  ∂z ∂ζ   ∂z [J ]： ヤコビのマトリクス（Jacobi matrix
Jacobian）

FEM3D34

自然座標系における偏微分（3/4）• Niの定義より簡単に求められる
∂y
∂  8
∂x
∂  8
 8 ∂N i 8 ∂N i
J11 =xi , J12 =yi ,==  N i yi  =   N i xi  = ∂ξ ∂ξ  i =1∂ξ ∂ξ  i =1
 i =1 ∂ξ i =1 ∂ξ∂z∂  8 8 ∂N i
J13 =zi=
  N i zi  = ∂ξ ∂ξ  i =1
 i =1 ∂ξ
∂x
∂  8
∂y
∂  8
 8 ∂N i 8 ∂N i
J 21 =
xi , J 22 =yi ,==
  N i xi  =   N i yi  = ∂η ∂η  i =1∂η ∂η  i =1
 i =1 ∂η i =1 ∂η∂z∂  8 8 ∂N i
J 23 =zi=  N i zi  = ∂η ∂η  i =1
 i =1 ∂η
∂x
∂  8
∂y
∂  8
 8 ∂N i 8 ∂N i
J 31 =
xi , J 32 =yi ,==
  N i xi  =   N i yi  = ∂ζ ∂ζ  i =1∂ζ ∂ζ  i =1
 i =1 ∂ζ i =1 ∂ζ∂z∂  8 8 ∂N i
J 33 ==zi  N i zi  = ∂ζ ∂ζ  i =1
 i =1 ∂ζ

FEM3D

35

自然座標系における偏微分（4/4）• 従って下記のように偏微分を計算できる– ヤコビアン（3×3行列）の逆行列を求める
 ∂N i   ∂x
 ∂x   ∂ξ ∂N   ∂x
 i = ∂y   ∂η ∂N i   ∂x ∂z   ∂ζ

∂y
∂ξ
∂y
∂η
∂y
∂ζ

∂z 
∂ξ 
∂z 
∂η 
∂z 

∂ζ 

−1

 ∂N i  ∂N i  ∂ξ  ∂ξ 
 ∂N i −1  ∂N i 
[
]
=
J




η
∂
∂
η




 ∂N i  ∂N i  ∂ζ  ∂ζ 

FEM3D36

要素単位での積分[k ] (i, j = 1…8)
ij(− 1,+1,+1) 8

7 (+ 1,+1,+1)(− 1,−1,+1) 5

6 (+ 1,−1,+1)(− 1,+1,−1)
4
1(ξ ,η , ζ ) = (− 1,−1,−1)kij = −  {λ ⋅ N i , x ⋅ N j , x + λ ⋅ N i , y ⋅ N j , y + λ ⋅ N i , z ⋅ N j , z }dVV ∂N i ∂N j
∂N i ∂N j∂N i ∂N j = −  λ
+λ
+λ
dV
∂x ∂x
∂y ∂y
∂z ∂z V 

2

3 (+ 1,+1,−1)(+ 1,−1,−1)

FEM3D37

自然座標系での積分[k ] (i, j = 1…8)
ij(− 1,+1,+1) 8(− 1,−1,+1) 5

7 (+ 1,+1,+1)
6 (+ 1,−1,+1)(− 1,+1,−1)

3 (+ 1,+1,−1)

4
1(ξ ,η , ζ ) = (− 1,−1,−1)2(+ 1,−1,−1) ∂N i ∂N j
∂N i ∂N j∂N i ∂N j −  λ
+λ
+λ
dV =
∂x ∂x
∂y ∂y
∂z ∂z V  ∂N i ∂N j
∂N i ∂N j∂N i ∂N j −    λ+λ+λ
dxdydz =∂y ∂y∂z ∂z  ∂x ∂x+1 +1 +1
 ∂N i ∂N j
∂N i ∂N j∂N i ∂N j −    λ+λ+λ
 det J dξdηdζ
∂x ∂x
∂y ∂y
∂z ∂z −1 −1 −1 

FEM3D38

ガウスの積分公式[k ] (i, j = 1…8)(− 1,+1,+1) 8

ij(− 1,−1,+1) 5

7 (+ 1,+1,+1)
6 (+ 1,−1,+1)(− 1,+1,−1)
4
1(ξ ,η , ζ ) = (− 1,−1,−1)

2

3 (+ 1,+1,−1)(+ 1,−1,−1) ∂N i ∂N j
∂N i ∂N j∂N i ∂N j −    λ+λ+λ
 det J dξdηdζ
∂x ∂x
∂y ∂y
∂z ∂z −1 −1 −1 
+1 +1 +1

+1 +1 +1

  f (ξ ,η , ζ ) dξdηdζ

I=

−1 −1 −1
L

M

N

i =1

j =1

k =1

[

=    Wi ⋅ W j ⋅ Wk ⋅ f (ξ i ,η j , ζ k )]FEM3D39

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

FEM3D

40

要素⇒全体マトリクス重ね合わせ
5

6

3[k ( 2 ) ]{φ ( 2) } = { f ( 2 ) }12

k11( 2)
 ( 2)
k 21
k31( 2)
 ( 2)
k 414

3[k (1) ]{φ (1) } = { f (1) }2

k11(1) (1)k 21
k31(1) (1)k 41

4
1

2
4

3
1

2
1

2
1

k12( 2)( 2)
k 22

k13( 2)( 2)
k 23

k32( 2)k33( 2)

( 2)
k 42

( 2)
k 43

k12(1)(1)
k 22

k13(1)(1)
k 23

k32(1)

k33(1)

(1)
k 42

(1)
k 43

k14( 2)  φ1( 2 )   f1( 2)  ( 2) ( 2)   ( 2) 
k 24
φ


f 
 2
=  2( 2 ) ( 2)   ( 2) 
k34 φ3   f 3 ( 2)( 2)   ( 2) k 44
 φ4   f 4 

k14(1)  φ1(1)   f1(1)  (1) (1)   (1) 
k 24
φ


f 
2

=  2(1) (1)   (1) 
k34 φ3   f 3 (1)(1)   (1) k 44
 φ4   f 4 

[ K ]{Φ} = {F }
 D1
X

X

X




X

X

X

D2
X

X
D3

X
X

X

X

X
X
X

D4
X
X

X
D5
X

  Φ1   B1  Φ   B 
 2   2 
X  Φ 3   B3   =  X  Φ 4   B4 
X  Φ 5   B5    
D6  Φ 6   B6 FEM3D41• 三次元要素の定式化• 三次元弾性力学方程式– ガラーキン法– 要素マトリクス生成– 演習• プログラムの実行• データ構造• プログラムの構成

FEM3D

42

演習• ガウスの積分公式を使用して以下の四角形の面積を求めよ（プログラムを作って，計算機で計算してください）y3
4

1:(1.0, 1.0)
2:(4.0, 2.0)
3:(3.0, 5.0)
4:(2.0, 4.0)

V
2
1
x
+1 +1

I =  dV = V det J dξdζ−1 −1

FEM3D

43

ヒント（1/2）• 座標値によってヤコビアン（ヤコビの行列）を計算• ガウスの積分公式（n=2）に代入する。
+1 +1

I=

  f (ξ ,η ) dξdη =   [W ⋅W ⋅ f (ξ ,η )]
m

n

i

−1 −1

i =1

j =1

implicit REAL*8 (A-H,O-Z)
real*8 W(2)real*8 POI(2)
W(1)= 1.0d0
W(2)= 1.0d0
POI(1)= -0.5773502692d0
POI(2)= +0.5773502692d0
SUM= 0.d0
do jp= 1, 2
do ip= 1, 2
FC = F(POI(ip),POI(jp))
SUM= SUM + W (ip)*W (jp)*FC
enddo
enddo

j

i

j

FEM3D

44

ヒント（2/2） ∂x ∂ξ[J ] =  ∂x

 ∂η

∂y 
∂x ∂y ∂y ∂x
∂ξ 
⋅
−
⋅
, det J =∂y ∂ξ ∂η ∂ξ ∂η
∂η 

∂x
∂  4
 4 ∂N ixi ,=  N i xi  = ∂ξ ∂ξ  i =1
 i =1 ∂ξ∂y∂  4 4 ∂N iyi ,=  N i yi  = ∂ξ ∂ξ  i =1
 i =1 ∂ξ∂x∂  4 4 ∂N ixi ,=  N i xi  = ∂η ∂η  i =1
 i =1 ∂η∂y∂  4 4 ∂N iyi=  N i yi  = ∂η ∂η  i =1
 i =1 ∂η1
1
N1 (ξ ,η ) = (1 − ξ )(1 − η ), N 2 (ξ ,η ) = (1 + ξ )(1 − η ),441
1
N 3 (ξ ,η ) = (1 + ξ )(1 + η ), N 4 (ξ ,η ) = (1 − ξ )(1 + η )4
4

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

FEM3D46

対象とする問題：三次元定常熱伝導∂  ∂T  ∂  ∂T  ∂  ∂T  ɺ +  λ + Q ( x, y , z ) = 0
λ +  λ∂x  ∂x  ∂y  ∂y  ∂z  ∂z 

• 定常熱伝導＋発熱• 一様な熱伝導率 λ
• 直方体

Z

T=0@Z=zmax– 一辺長さ1の立方体（六面体）要素– 各方向にNX・NY・NZ個
NZ• 境界条件– T=0@Z=zmax

Y
NX
NY
X

• 体積当たり発熱量は位置（メッシュの中心の座標 xc,yc）に依存– Qɺ ( x, y, z ) = QVOL xC + yC

FEM3D47

インストール（Cygwinでは*.exe）
インストール>$ cd <$P-TOP>/fem3d/src>$ make>$ ls ../run/sol
sol
メッシュジェネレータインストール>$ cd <$P-TOP>/fem3d/run>$ cc –O3 mgcube.c –o mgcube

FEM3D48

計算の流れ
メッシュ生成⇒計算，ファイル名称固定
mgcubeメッシュジェネレータ

cube.0sol

メッシュファイル

FEM本体

test.inp
可視化用出力INPUT.DAT
制御データ

FEM3D49

メッシュ生成（Cygwinではmgcube.exe）>$ cd <$P-TOP>/fem3d/run>$ ./mgcube

Z

T=0@Z=zmax

NX, NY, NZ←

20,20,20←>$ ls cube.0
NZ

cube.0
Y
NX
NY
X

各辺長さを
訊いてくる
このように
入れてみる
生成を確認

FEM3D

50

制御ファイル：INPUT.DAT
INPUT.DAT
cube.02000
1.0 1.0
1.0e-08• fname：• ITER：• COND：• QVOL：• RESID：

fname
ITER
COND, QVOL
RESID

メッシュファイル名
反復回数上限熱伝導率
体積当たり発熱量係数
反復法の収束判定値∂  ∂T  ∂  ∂T  ∂  ∂T  ɺ +  λλ +  λ + Q ( x, y , z ) = 0∂x  ∂x  ∂y  ∂y  ∂z  ∂z Qɺ ( x, y, z ) = QVOL xC + yC

FEM3D51

実 行（Cygwinではsol.exe）>$ cd <$P-TOP>/fem3d/run>$ ./sol>$ ls test.inp
test.inp

生成を確認

FEM3D

52

ParaView• http://www.paraview.org/• ファイルを開く• 図の表示• イメージファイルの保存• http://nkl.cc.u-tokyo.ac.jp/class/ParaView.pdf

FEM3D53

UCD Format (1/3)
Unstructured Cell Data
要素の種類

キーワード

点

pt

線

line

三角形

tri

四角形

quad

四面体

tet

角錐

pyr

三角柱

prism

六面体

hex

二次要素
線2

line2

三角形2

tri2

四角形2

quad2

四面体2

tet2

角錐2

pyr2

三角柱2

prism2六面体2

hex2

FEM3D54

UCD Format (2/3)• Originally for AVS, microAVS• Extension of the UCD file is “inp”• There are two types of formats. Only old type canbe read by ParaView.

FEM3D55

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

・
・
・

・
・
・

(節点のデータ成分数) (成分1の構成数) (成分2の構成数) ・・・(各成分の構成数)(節点データ成分1のラベル)，(単位)(節点データ成分2のラベル)，(単位)・・・(各節点データ成分のラベル)，(単位)(節点番号1) (節点データ1) (節点データ2) ・・・・・(節点番号2) (節点データ1) (節点データ2) ・・・・・
・・・

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

FEM3D57

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

cube.0：節点データ（NX=NY=NZ=4）
125
1
2
3
4
5
6
7
8
9

0.00
1.00
2.00
3.00
4.00
0.00
1.00
2.00
3.00

0.00
0.00
0.00
0.00
0.00
1.00
1.00
1.00
1.00

0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00

121
122
123
124
125

0.00
1.00
2.00
3.00
4.00

4.00
4.00
4.00
4.00
4.00

4.00
4.00
4.00
4.00
4.00

節点ID

X座標

Y座標

Z座標

=5*5*5（節点数）Z

T=0@Z=zmax

NZ

...
Y
NX
NY
X

movie

FEM3D59

cube.0：要素データ（1/2）
64=4*4*4（要素数）361 361 361 361 361 361 361 361 361 361
361 361 361 361 361 361 361 361 361 361
361 361 361 361 361 361 361 361 361 361
361 361 361 361 361 361 361 361 361 361
361 361 361 361 361 361 361 361 361 361
361 361 361 361 361 361 361 361 361 361
361 361 361 361

Z

T=0@Z=zmax

NZ

要素タイプ：361
三次元，六面体，一次

Y
NX
NY
X

FEM3D60

cube.0：要素データ（2/2）

...

1
2
3
4
5
6
7
8
9
10
11
12
13

1
1
1
1
1
1
1
1
1
1
1
1
1

1
2
3
4
6
7
8
9
11
12
13
14
16

2
3
4
5
7
8
9
10
12
13
14
15
17

53
54
55
56
57
58
59
60
61
62
63
64

1
1
1
1
1
1
1
1
1
1
1
1

81
82
83
84
86
87
88
89
91
92
93
94

82 87
83 88
84 89
85 90
87 92
88 93
89 94
90 95
92 97
93 98
94 99
95 100

要素ID 材料ID

7
8
9
10
12
13
14
15
17
18
19
20
22

6
7
8
9
11
12
13
14
16
17
18
19
21

26
27
28
29
31
32
33
34
36
37
38
39
41

27
28
29
30
32
33
34
35
37
38
39
40
42

32
33
34
35
37
38
39
40
42
43
44
45
47

86 106 107 112 111
87 107 108 113 112
88 108 109 114 113
89 109 110 115 114
91 111 112 117 116
92 112 113 118 117
93 113 114 119 118
94 114 115 120 119
96 116 117 122 121
97 117 118 123 122
98 118 119 124 123
99 119 120 125 124

8節点ID

Z

31
32
33
34
36
37
38
39
41
42
43
44
46

T=0@Z=zmax

NZ

Y
NX
NY
X

8
5

7
6

4
1

3
2

FEM3D61

NX=NY=NZ=4, NXP1=NYP1=NZP1=5
ICELTOT= 64, INODTOT= 125, IBNODTOT= 25
k=2

k=1
24

25
20

1
1

8

2
2

21

9

3
3

10 31

4
5

26

k=5

34

18

19

20

27

28

106

107

108

109

50

51

52

49
102

103

30

51

104

34

35

36

53

64
95

87

110 81

82

83

84

50

51

52

49
76

59

100

94

53

105

58

52

99

86

77

70

57

33

120

112

53

56

64

119
111

35

29

48

62

37

125

124

101

33

75

69
61

32

17

4

45

37

36

7

6

74

32

44

12

5

50

49

16

19
11

k=3

78

79

85

80

k=4

54

60

55

FEM3D62

cube.0：節点グループデータ
4
25 50

グループ総数
各グループ節点数（累積）75 100Xmin
1 6 11 16 21
51 56 61 66 71
101 106 111 116 121

26
76

31
81

36
86

41
91

46
96

Ymin
1 2 3 4 551 52 53 54 55
101 102 103 104 105

26
76

27
77

28
78

29
79

30
80

Zmin
1 2
11 12
21 22

T=0@Z=zmax

NZ

Y

3
13
23

4
14
24

5
15
25

6
16

7
17

8
18

9
19

10
20

Zmax
101 102 103 104 105 106 107 108 109 110
111 112 113 114 115 116 117 118 119 120
121 122 123 124 125
（以下

Z

使用せず）

NX
NY
X

FEM3D63

NX=NY=NZ=4, NXP1=NYP1=NZP1=5
ICELTOT= 64, INODTOT= 125, IBNODTOT= 25
k=2

k=1
24

25
20

1
1

8

2
2

21

9

3
3

10 31

4
5

26

k=5

34

18

19

20

27

28

106

108

109

50

51

52

49
102

103

30

51

104

34

35

36

53

64

54

60

55

k=4

95

87

110 81

82

83

84

50

51

52

49
76

59

100

94

53

105

58

52

99

86

107

70

57

33

120

112

53

56

64

119
111

35

29

48

62

37

125

124

101

33

75

69
61

32

17

4

45

37

36

7

6

74

32

44

12

5

50

49

16

19
11

k=3

77

78

79

85

80

Xmin: i=1
Ymin: j=1
Zmin: k=1
Zmax:k=5

FEM3D

64

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

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

FEM3D66

有限要素法の処理：プログラム• 初期化– 制御変数読み込み– 座標読み込み⇒要素生成（N:節点数，NE：要素数）– 配列初期化（全体マトリクス，要素マトリクス）– 要素⇒全体マトリクスマッピング（Index，Item）• マトリクス生成– 要素単位の処理（do icel= 1, NE）• 要素マトリクス計算• 全体マトリクスへの重ね合わせ– 境界条件の処理• 連立一次方程式– 共役勾配法（CG）

FEM3D67

test1

input_cntl

メインプログラム制御データ入力

input_grid

find_node

メッシュファイル入力

節点探索

mat_con0mSORT

行列コネクティビティ生成

ソート

mat_con1
行列コネクティビティ生成

mat_ass_main

jacobi係数行列生成ヤコビアン計算

mat_ass_bc
境界条件処理三次元熱伝導解析コードheat3Dの構成

solve11cg

線形ソルバー制御

CG法計算

output_ucd
可視化処理FEM3D

68

全体処理
/**

program heat3D**/#include <stdio.h>#include <stdlib.h>
FILE* fp_log;#define GLOBAL_VALUE_DEFINE#include "pfem_util.h"//#include "solver11.h"extern void INPUT_CNTL();
extern void INPUT_GRID();
extern void MAT_CON0();
extern void MAT_CON1();
extern void MAT_ASS_MAIN();
extern void MAT_ASS_BC();
extern void SOLVE11();
extern void OUTPUT_UCD();
int main(){
INPUT_CNTL();INPUT_GRID();
MAT_CON0();
MAT_CON1();
MAT_ASS_MAIN();
MAT_ASS_BC() ;SOLVE11();}OUTPUT_UCD();FEM3D

69

Global変数表：pfem_util.h（1/3）
変数名

種別

サイズ

内 容

I/O
I

メッシュファイル名

I

I

節点数

ICELTOTII

要素数NODGRPtot

I

I

節点グループ数

XYZ

R[N][3]I

節点座標

ICELNODI[ICELTOT][8]I

要素コネクティビティ

NODGRP_INDEXI[NODGRPtot+1]I

各節点グループに含まれる節点数（累積）NODGRP_ITEMI[NODGRP_INDEX[NODG
RPTOT+1]]

I

節点グループに含まれる節点

NODGRP_NAME

C80[NODGRP_INDEX[NODG
RPTOT+1]]

I

節点グループ名

NLU

IO

各節点非対角成分数

NPLU

I

O

非対角成分総数

D

R

[N]

O

全体行列：対角ブロック

B,X

R

[N]

O

右辺ベクトル，未知数ベクトル

fname

C

N,NP

[80]

FEM3D

70

Global変数表：pfem_util.h（2/3）
変数名

種別

サイズ

I/O

内 容

AMAT

R

[N]

O

全体行列：非零非対角成分

indexLU

I

[N+1]

O

全体行列：非零非対角成分数

itemLUI[NPLU]

O

全体行列：非零非対角成分（列番号）

INLU

I

[N]

O

各節点の非零非対角成分数

IALU

I

[N][NLU]

O

各節点の非零非対角成分数（列番号）

IWKX

I

[N][2]

O

ワーク用配列ITER, ITERactualI

I

反復回数の上限，実際の反復回数

RESID

R

I

打ち切り誤差（1.e-8に設定）

pfemIarray

I

[100]

O

諸定数（整数）

pfemRarray

R

[100]

O

諸定数（実数）FEM3D

71

Global変数表：pfem_util.h（3/3）
変数名

種別

O8th

R

PNQ, PNE, PNT

R

POS, WEI

サイズ

内 容

I/O
I

=0.125[2][2][8]O

各ガウス積分点における

R

[2]

O

各ガウス積分点の座標，重み係数

NCOL1, NCOL2

I

[100]

O

ソート用ワーク配列

SHAPE

R

[2][2][2][8]O

各ガウス積分点における形状関数 Ni (i=1~8)

PNX, PNY, PNZ

R[2][2][2][8]O

各ガウス積分点における

DETJ

R

[2][2][2]O

各ガウス積分点におけるヤコビアン行列式

COND, QVOL

R

I

熱伝導率，体積当たり発熱量係数∂N i ∂N i ∂N i
(i = 1 ~ 8),,∂ξ ∂η ∂ζ∂N i ∂N i ∂N i
(i = 1 ~ 8),,∂x ∂y ∂z

FEM3D

72

制御ファイル入力：INPUT_CNTL/**** INPUT_CNTL**/#include <stdio.h>#include <stdlib.h>#include "pfem_util.h"/** **/
void INPUT_CNTL(){
FILE *fp;if( (fp=fopen("INPUT.DAT","r")) == NULL){
fprintf(stdout,"input file cannot be opened!¥n");exit(1);}
fscanf(fp,"%s",fname);
INPUT.DAT
fscanf(fp,"%d",&ITER);
fscanf(fp, "%lf %lf", &COND, &QVOL);
fscanf(fp, "%lf",&RESID);
cube.0fclose(fp);2000}

pfemRarray[0]= RESID;
pfemIarray[0]= ITER;

1.0 1.0
1.0e-08

fname
ITER
COND, QVOL
RESID

73

FEM3D

メッシュ入力：INPUT_GRID（1/3）#include <stdio.h>#include <stdlib.h>#include "pfem_util.h"#include "allocate.h"void INPUT_GRID(){
FILE *fp;
int i,j,k,ii,kk,nn,icel,iS,iE;
int NTYPE,IMAT;

/**
**/

if( (fp=fopen(fname,"r")) == NULL){
fprintf(stdout,"input file cannot be opened!¥n");exit(1);}NODE
fscanf(fp,"%d",&N);
NP=N;
XYZ=(KREAL**)allocate_matrix(sizeof(KREAL),N,3);
for(i=0;i<N;i++){for(j=0;j<3;j++){
XYZ[i][j]=0.0;}}
for(i=0;i<N;i++){fscanf(fp,"%d %lf %lf %lf",&ii,&XYZ[i][0],&XYZ[i][1],&XYZ[i][2]);}FEM3D74

cube.0：節点データ（NX=NY=NZ=4）
125
1
2
3
4
5
6
7
8
9

0.00
1.00
2.00
3.00
4.00
0.00
1.00
2.00
3.00

0.00
0.00
0.00
0.00
0.00
1.00
1.00
1.00
1.00

0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00

0.00
1.00
2.00
3.00
4.00

4.00
4.00
4.00
4.00
4.00

4.00
4.00
4.00
4.00
4.00

= N

Z

T=0@Z=zmax

NZ

...
121
122
123
124
125

XYZ[i][3]

Y
NX
NY
X

FEM3D

allocate, deallocate (1/2)
allocateをFORTRAN並みに簡単にやるための関数#include <stdio.h>#include <stdlib.h>
void* allocate_vector(int size,int m){void *a;if ( ( a=(void * )malloc( m * size ) ) == NULL ) {fprintf(stdout,"Error:Memory does not enough! in vector ¥n");exit(1);}
return a;}void deallocate_vector(void *a){free( a );}INDEX=(KINT* )allocate_vector(sizeof(KINT),NGtot+1); INDEX[NGtot]
NAME =(CHAR80*)allocate_vector(sizeof(CHAR80),NGtot); NAME[NGtot]
WW=(KREAL**) allocate_matrix(sizeof(KREAL),4,N);
WW[4][N]75FEM3D

allocate, deallocate (2/2)
allocateをFORTRAN並みに簡単にやるための関数

void** allocate_matrix(int size,int m,int n){void **aa;
int i;
if ( ( aa=(void ** )malloc( m * sizeof(void*) ) ) == NULL ) {fprintf(stdout,"Error:Memory does not enough! aa in matrix ¥n");exit(1);}
if ( ( aa[0]=(void * )malloc( m * n * size ) ) == NULL ) {fprintf(stdout,"Error:Memory does not enough! in matrix ¥n");exit(1);}for(i=1;i<m;i++) aa[i]=(char*)aa[i-1]+size*n;
return aa;}
void deallocate_matrix(void **aa){
free( aa );}INDEX=(KINT* )allocate_vector(sizeof(KINT),NGtot+1); INDEX[NGtot]
NAME =(CHAR80*)allocate_vector(sizeof(CHAR80),NGtot); NAME[NGtot]
WW=(KREAL**) allocate_matrix(sizeof(KREAL),4,N);
WW[4][N]76FEM3D77

メッシュ入力：INPUT_GRID（2/3）/****/

ELEMENTICELNOD[i][j]の中身としては
fscanf(fp,"%d",&ICELTOT);「1」から始まる通し節点番号がICELNOD=(KINT**)allocate_matrix(sizeof(KINT),ICELTOT,8); そのまま読み込まれている。
for(i=0;i<ICELTOT;i++) fscanf(fp,"%d",&NTYPE);
要素番号は「0」から番号付け。
for(icel=0;icel<ICELTOT;icel++){fscanf(fp,"%d %d %d %d %d %d %d %d %d %d",&ii,&IMAT,&ICELNOD[icel][0],&ICELNOD[icel][1],&ICELNOD[icel][2],&ICELNOD[icel][3],&ICELNOD[icel][4],&ICELNOD[icel][5],&ICELNOD[icel][6],&ICELNOD[icel][7]);}FEM3D78

cube.0：要素データ（1/2）
64= ICELTOTZ

361 361 361 361 361 361 361 361 361 361
361 361 361 361 361 361 361 361 361 361
361 361 361 361 361 361 361 361 361 361
361 361 361 361 361 361 361 361 361 361
361 361 361 361 361 361 361 361 361 361
361 361 361 361 361 361 361 361 361 361
361 361 361 361

T=0@Z=zmax

NZ

要素タイプ：361
三次元，六面体，一次

Y
NX
NY
X

FEM3D79

cube.0：要素データ（2/2）

...

1
2
3
4
5
6
7
8
9
10
11
12
13

1
1
1
1
1
1
1
1
1
1
1
1
1

1
2
3
4
6
7
8
9
11
12
13
14
16

2
3
4
5
7
8
9
10
12
13
14
15
17

53
54
55
56
57
58
59
60
61
62
63
64

1
1
1
1
1
1
1
1
1
1
1
1

81
82
83
84
86
87
88
89
91
92
93
94

82 87
83 88
84 89
85 90
87 92
88 93
89 94
90 95
92 97
93 98
94 99
95 100

iMAT

7
8
9
10
12
13
14
15
17
18
19
20
22

6
7
8
9
11
12
13
14
16
17
18
19
21

26
27
28
29
31
32
33
34
36
37
38
39
41

27
28
29
30
32
33
34
35
37
38
39
40
42

32
33
34
35
37
38
39
40
42
43
44
45
47

86 106 107 112 111
87 107 108 113 112
88 108 109 114 113
89 109 110 115 114
91 111 112 117 116
92 112 113 118 117
93 113 114 119 118
94 114 115 120 119
96 116 117 122 121
97 117 118 123 122
98 118 119 124 123
99 119 120 125 124

ICELNOD[icel][8]

Z

31
32
33
34
36
37
38
39
41
42
43
44
46

T=0@Z=zmax

NZ

Y
NX
NY
X

8
5

7
6

4
1

3
2

FEM3D80

メッシュ入力：INPUT_GRID（3/3）/****/

NODE grp. info.
fscanf(fp,"%d",&NODGRPtot);
NODGRP_INDEX=(KINT* )allocate_vector(sizeof(KINT),NODGRPtot+1);
NODGRP_NAME =(CHAR80*)allocate_vector(sizeof(CHAR80),NODGRPtot);for(i=0;i<NODGRPtot+1;i++) NODGRP_INDEX[i]=0;
for(i=0;i<NODGRPtot;i++) fscanf(fp,"%d",&NODGRP_INDEX[i+1]);
nn=NODGRP_INDEX[NODGRPtot];
NODGRP_ITEM=(KINT*)allocate_vector(sizeof(KINT),nn);for(k=0;k<NODGRPtot;k++){iS= NODGRP_INDEX[k];iE= NODGRP_INDEX[k+1];
fscanf(fp,"%s",NODGRP_NAME[k].name);
nn= iE - iS;
if( nn != 0 ){for(kk=iS;kk<iE;kk++) fscanf(fp,"%d",&NODGRP_ITEM[kk]);}}}

fclose(fp);

節点グループの中身も「1」から始まる通し節点番号がそのまま読み込まれている。

FEM3D81

cube.0：節点グループデータ
425 50

NODGRPtot
NODGTP_INDEX[1-4]

75 100Xmin
1 6 11 16 21 26 31 36 41 46
51 56 61 66 71 76 81 86 91 96
101 106 111 116 121

NODGRP_NAME[1]NODGRP_ITEM[0-24]

Ymin
1 2 3 45 26 27 28 29 30
51 52 53 54 55 76 77 78 79 80
101 102 103 104 105

NODGRP_NAME[2]NODGRP_ITEM[25-49]

Zmin
1 2
11 12
21 22

NODGRP_NAME[3]NODGRP_ITEM[50-74]

3
13
23

4
14
24

5 6 7 8 9 1015 16 17 18 19 20
25

Zmax
101 102 103 104 105 106 107 108 109 110
111 112 113 114 115 116 117 118 119 120
121 122 123 124 125

NODGRP_NAME[4]NODGRP_ITEM[75-99]

FEM3D82

test1

input_cntl

メインプログラム制御データ入力

input_grid

find_node

メッシュファイル入力

節点探索

mat_con0mSORT

行列コネクティビティ生成

ソート

mat_con1
行列コネクティビティ生成

mat_ass_main

jacobi係数行列生成ヤコビアン計算

mat_ass_bc
境界条件処理三次元熱伝導解析コードheat3Dの構成

solve11cg

線形ソルバー制御

CG法計算

output_ucd
可視化処理FEM3D

83

Global変数表：pfem_util.h（1/3）
変数名

種別

サイズ

内 容

I/O
I

メッシュファイル名

I

I

節点数

ICELTOTII

要素数NODGRPtot

I

I

節点グループ数

XYZ

R[N][3]I

節点座標

ICELNODI[ICELTOT][8]I

要素コネクティビティ

NODGRP_INDEXI[NODGRPtot+1]I

各節点グループに含まれる節点数（累積）NODGRP_ITEMI[NODGRP_INDEX[NODG
RPTOT+1]]

I

節点グループに含まれる節点

NODGRP_NAME

C80[NODGRP_INDEX[NODG
RPTOT+1]]

I

節点グループ名

NLU

IO

各節点非対角成分数

NPLU

I

O

非対角成分総数

D

R

[N]

O

全体行列：対角ブロック

B,X

R

[N]

O

右辺ベクトル，未知数ベクトル

fname

C

N,NP

[80]

FEM3D

84

Global変数表：pfem_util.h（2/3）
変数名

種別

サイズ

I/O

内 容

AMAT

R

[N]

O

全体行列：非零非対角成分

indexLU

I

[N+1]

O

全体行列：非零非対角成分数

itemLUI[NPLU]

O

全体行列：非零非対角成分（列番号）

INLU

I

[N]

O

各節点の非零非対角成分数

IALU

I

[N][NLU]

O

各節点の非零非対角成分数（列番号）

IWKX

I

[N][2]

O

ワーク用配列ITER, ITERactualI

I

反復回数の上限，実際の反復回数

RESID

R

I

打ち切り誤差（1.e-8に設定）

pfemIarray

I

[100]

O

諸定数（整数）

pfemRarray

R

[100]

O

諸定数（実数）FEM3D

85

Global変数表：pfem_util.h（3/3）
変数名

種別

O8th

R

PNQ, PNE, PNT

R

POS, WEI

サイズ

内 容

I/O
I

=0.125[2][2][8]O

各ガウス積分点における

R

[2]

O

各ガウス積分点の座標，重み係数

NCOL1, NCOL2

I

[100]

O

ソート用ワーク配列

SHAPE

R

[2][2][2][8]O

各ガウス積分点における形状関数 Ni (i=1~8)

PNX, PNY, PNZ

R[2][2][2][8]O

各ガウス積分点における

DETJ

R

[2][2][2]O

各ガウス積分点におけるヤコビアン行列式

COND, QVOL

R

I

熱伝導率，体積当たり発熱量係数∂N i ∂N i ∂N i
(i = 1 ~ 8),,∂ξ ∂η ∂ζ∂N i ∂N i ∂N i
(i = 1 ~ 8),,∂x ∂y ∂z∂  ∂T  ∂  ∂T  ∂  ∂T  ɺ +  λλ +  λ + Q ( x, y , z ) = 0∂x  ∂x  ∂y  ∂y  ∂z  ∂z Qɺ ( x, y, z ) = QVOL xC + yC

FEM3D86

マトリクス生成まで• 一次元のときは，index，itemに関連した情報を簡単に作ることができた– 非ゼロ非対角成分の数は2– 番号が自分に対して：+1と-1• 三次元の場合はもっと複雑– 非ゼロ非対角成分の数は7～26（現在の形状）– 実際はもっと複雑– 前以て，非ゼロ非対角成分数はわからないmovie

FEM3D87

マトリクス生成まで• 一次元のときは，index，itemに関連した情報を簡単に作ることができた– 非ゼロ非対角成分の数は2– 番号が自分に対して：+1と-1• 三次元の場合はもっと複雑– 非ゼロ非対角成分の数は7～26（現在の形状）– 実際はもっと複雑– 前以て，非ゼロ非対角成分数はわからない• INLU[N]，IALU[N][NLU] を使って非ゼロ非対角成分数を予備的に勘定する

FEM3D

88

全体処理
/**

program heat3D**/#include <stdio.h>#include <stdlib.h>
FILE* fp_log;#define GLOBAL_VALUE_DEFINE#include "pfem_util.h"//#include "solver11.h"extern void INPUT_CNTL();
extern void INPUT_GRID();
extern void MAT_CON0();
extern void MAT_CON1();
extern void MAT_ASS_MAIN();
extern void MAT_ASS_BC();
extern void SOLVE11();
extern void OUTPUT_UCD();
int main()
{

13

14

7
9

8
10

4
5

16

9
11

5
6

1
1

15

12

6
7

2
2

8

3
3

4

INPUT_CNTL();INPUT_GRID();

MAT_CON0: INU, IALU生成
MAT_CON1: index, item生成

MAT_CON0();
MAT_CON1();
MAT_ASS_MAIN();
MAT_ASS_BC() ;SOLVE11();}OUTPUT_UCD();

とりあえず1から始まる節点番号を記憶

FEM3D89

MAT_CON0：全体構成(− 1,+1,+1) 8

do icel= 1, ICELTOT
8節点相互の関係から，
INLU，IALUを生成（FIND_NODE）(− 1,−1,+1) 5

7 (+ 1,+1,+1)
6 (+ 1,−1,+1)

enddo(− 1,+1,−1)
4
1

13

14

7
9

8
10

4
5

11

6

12

6
7

2
2

16

9

5

1
1

15

8

3
3

4

(ξ ,η , ζ ) = (− 1,−1,−1)

2

3 (+ 1,+1,−1)(+ 1,−1,−1)

FEM3D90

行列コネクティビティ生成：
MAT_CON0（1/4）/**** MAT_CON0**/#include <stdio.h>#include "pfem_util.h"#include "allocate.h“
extern FILE *fp_log;/*** external functions ***/
extern void mSORT(int*, int*, int);/*** static functuons ***/
static void FIND_TS_NODE (int,int);
void MAT_CON0(){
int i,j,k,icel,in;
int in1,in2,in3,in4,in5,in6,in7,in8;
int NN;NLU= 26;
INLU=(KINT* )allocate_vector(sizeof(KINT),N);IALU=(KINT**)allocate_matrix(sizeof(KINT),N,NLU);
for(i=0;i<N;i++) INLU[i]=0;
for(i=0;i<N;i++) for(j=0;j<NLU;j++) IALU[i][j]=0;

NLU：
各節点における
非ゼロ非対角成分の最大数（接続する節点数）
今の問題の場合はわかっているので，
このようにできる不明の場合の実装：⇒レポート課題

FEM3D91

行列コネクティビティ生成：
MAT_CON0（1/4）/**** MAT_CON0**/#include <stdio.h>#include "pfem_util.h"#include "allocate.h“
extern FILE *fp_log;/*** external functions ***/
extern void mSORT(int*, int*, int);/*** static functuons ***/
static void FIND_TS_NODE (int,int);

変数名

サイズ

内

容

INLU

[N]

各節点の非零非対角成
分数

IALU

[N]{[NLU]

各節点の非零非対角成分（列番号）

void MAT_CON0(){
int i,j,k,icel,in;
int in1,in2,in3,in4,in5,in6,in7,in8;
int NN;NLU= 26;
INLU=(KINT* )allocate_vector(sizeof(KINT),N);IALU=(KINT**)allocate_matrix(sizeof(KINT),N,NLU);
for(i=0;i<N;i++) INLU[i]=0;
for(i=0;i<N;i++) for(j=0;j<NLU;j++) IALU[i][j]=0;

FEM3D92

行列コネクティビティ生成：
MAT_CON0（2/4）：1から始まる番号
for( icel=0;icel< ICELTOT;icel++){in1=ICELNOD[icel][0];in2=ICELNOD[icel][1];in3=ICELNOD[icel][2];in4=ICELNOD[icel][3];in5=ICELNOD[icel][4];in6=ICELNOD[icel][5];in7=ICELNOD[icel][6];in8=ICELNOD[icel][7];FIND_TS_NODE (in1,in2);
FIND_TS_NODE (in1,in3);
FIND_TS_NODE (in1,in4);
FIND_TS_NODE (in1,in5);
FIND_TS_NODE (in1,in6);
FIND_TS_NODE (in1,in7);
FIND_TS_NODE (in1,in8);
FIND_TS_NODE (in2,in1);
FIND_TS_NODE (in2,in3);
FIND_TS_NODE (in2,in4);
FIND_TS_NODE (in2,in5);
FIND_TS_NODE (in2,in6);
FIND_TS_NODE (in2,in7);
FIND_TS_NODE (in2,in8);
FIND_TS_NODE (in3,in1);
FIND_TS_NODE (in3,in2);
FIND_TS_NODE (in3,in4);
FIND_TS_NODE (in3,in5);
FIND_TS_NODE (in3,in6);
FIND_TS_NODE (in3,in7);
FIND_TS_NODE (in3,in8);(− 1,+1,+1) 8(− 1,−1,+1) 5

7 (+ 1,+1,+1)
6 (+ 1,−1,+1)(− 1,+1,−1)
4
1(ξ ,η , ζ ) = (− 1,−1,−1)

2

3 (+ 1,+1,−1)(+ 1,−1,−1)

FEM3D93

節点探索：FIND_TS_NODEINLU,IAU探索：一次元ではこの部分は手動/****** FIND_TS_NODE**/
static void FIND_TS_NODE (int ip1,int ip2){
int kk,icou;
for(kk=1;kk<=INLU[ip1-1];kk++){if(ip2 == IALU[ip1-1][kk-1]) return;}icou=INLU[ip1-1]+1;
IALU[ip1-1][icou-1]=ip2;
INLU[ip1-1]=icou;

変数名

サイズ

return;
}

内

容

INLU

[N]

各節点の非零非対角成
分数

IALU

[N]{[NLU]

各節点の非零非対角成分（列番号）

FEM3D94

節点探索：FIND_TS_NODE
一次元ではこの部分は手動，要素#2/****** FIND_TS_NODE**/
static void FIND_TS_NODE (int ip1,int ip2){
int kk,icou;
for(kk=1;kk<=INLU[ip1-1];kk++){if(ip2 == IALU[ip1-1][kk-1]) return;}icou=INLU[ip1-1]+1;
IALU[ip1-1][icou-1]=ip2;
INLU[ip1-1]=icou;

ip2がIALU[ip1-1][kk]に含まれているかチェック

13

7
9

10

5

ip1: No.6 node
ip2: No.7 node

11

6

12

6
7

2
2

16

9

5

1
1

15

8

4

return;
}

14

8

3
3

4

FEM3D

95

節点探索：FIND_TS_NODE
一次元ではこの部分は手動，要素#2/****** FIND_TS_NODE**/
static void FIND_TS_NODE (int ip1,int ip2){
int kk,icou;
for(kk=1;kk<=INLU[ip1-1];kk++){if(ip2 == IALU[ip1-1][kk-1]) return;}icou=INLU[ip1-1]+1;
IALU[ip1-1][icou-1]=ip2;
INLU[ip1-1]=icou;

IALUに含まれていない
場合は，INLUに1を加えてIALUに格納
13

7
9

10

5

ip1: No.6 node
ip2: No.7 node

11

6

12

6
7

2
2

16

9

5

1
1

15

8

4

return;
}

14

8

3
3

4

FEM3D

96

節点探索：FIND_TS_NODE
一次元ではこの部分は手動，要素#5/****** FIND_TS_NODE**/
static void FIND_TS_NODE (int ip1,int ip2){
int kk,icou;
for(kk=1;kk<=INLU[ip1-1];kk++){if(ip2 == IALU[ip1-1][kk-1]) return;}icou=INLU[ip1-1]+1;
IALU[ip1-1][icou-1]=ip2;
INLU[ip1-1]=icou;

既にIALUに含まれている場合は，次のペアへ

13

7
9

10

5

ip1: No.6 node
ip2: No.7 node

11

6

12

6
7

2
2

16

9

5

1
1

15

8

4

return;
}

14

8

3
3

4

FEM3D97

行列コネクティビティ生成：
MAT_CON0（3/4）
FIND_TS_NODE (in4,in1);
FIND_TS_NODE (in4,in2);
FIND_TS_NODE (in4,in3);
FIND_TS_NODE (in4,in5);
FIND_TS_NODE (in4,in6);
FIND_TS_NODE (in4,in7);
FIND_TS_NODE (in4,in8);
FIND_TS_NODE (in5,in1);
FIND_TS_NODE (in5,in2);
FIND_TS_NODE (in5,in3);
FIND_TS_NODE (in5,in4);
FIND_TS_NODE (in5,in6);
FIND_TS_NODE (in5,in7);
FIND_TS_NODE (in5,in8);
FIND_TS_NODE (in6,in1);
FIND_TS_NODE (in6,in2);
FIND_TS_NODE (in6,in3);
FIND_TS_NODE (in6,in4);
FIND_TS_NODE (in6,in5);
FIND_TS_NODE (in6,in7);
FIND_TS_NODE (in6,in8);
FIND_TS_NODE (in7,in1);
FIND_TS_NODE (in7,in2);
FIND_TS_NODE (in7,in3);
FIND_TS_NODE (in7,in4);
FIND_TS_NODE (in7,in5);
FIND_TS_NODE (in7,in6);
FIND_TS_NODE (in7,in8);(− 1,+1,+1) 8(− 1,−1,+1) 5

7 (+ 1,+1,+1)
6 (+ 1,−1,+1)(− 1,+1,−1)
4
1(ξ ,η , ζ ) = (− 1,−1,−1)

2

3 (+ 1,+1,−1)(+ 1,−1,−1)

FEM3D98

行列コネクティビティ生成：
MAT_CON0（4/4）}

FIND_TS_NODE (in8,in1);
FIND_TS_NODE (in8,in2);
FIND_TS_NODE (in8,in3);
FIND_TS_NODE (in8,in4);
FIND_TS_NODE (in8,in5);
FIND_TS_NODE (in8,in6);
FIND_TS_NODE (in8,in7);

for(in=0;in<N;in++){
NN=INLU[in];for (k=0;k<NN;k++){
NCOL1[k]=IALU[in][k];}各節点において，
IALU[i][k]が小さい番号から大きい番号に並ぶようにソート（単純なバブルソート）
せいぜい100程度のものをソートする

mSORT(NCOL1,NCOL2,NN);}for(k=NN;k>0;k--){IALU[in][NN-k]= NCOL1[NCOL2[k-1]-1];}

FEM3D99

CRS形式への変換：MAT_CON1#include <stdio.h>#include "pfem_util.h"#include "allocate.h"extern FILE* fp_log;
void MAT_CON1(){
int i,k,kk;
indexLU=(KINT*)allocate_vector(sizeof(KINT),N+1);
for(i=0;i<N+1;i++) indexLU[i]=0;
for(i=0;i<N;i++){indexLU[i+1]=indexLU[i]+INLU[i];}Ciindex[i + 1] =  INLU[k ]k =0

NPLU=indexLU[N];itemLU=(KINT*)allocate_vector(sizeof(KINT),NPLU); index[0] = 0}

for(i=0;i<N;i++){for(k=0;k<INLU[i];k++){kk=k+indexLU[i];itemLU[kk]=IALU[i][k]-1;}}
deallocate_vector(INLU);
deallocate_vector(IALU);FORTRANiindex (i ) =  INLU(k )k =1

index (0) = 0

FEM3D

100

CRS形式への変換：MAT_CON1#include <stdio.h>#include "pfem_util.h"#include "allocate.h"extern FILE* fp_log;
void MAT_CON1(){
int i,k,kk;
indexLU=(KINT*)allocate_vector(sizeof(KINT),N+1);
for(i=0;i<N+1;i++) indexLU[i]=0;
for(i=0;i<N;i++){indexLU[i+1]=indexLU[i]+INLU[i];}NPLU=indexLU[N];itemLU=(KINT*)allocate_vector(sizeof(KINT),NPLU);}

for(i=0;i<N;i++){for(k=0;k<INLU[i];k++){kk=k+indexLU[i];itemLU[kk]=IALU[i][k]-1;}}
deallocate_vector(INLU);
deallocate_vector(IALU);NPLU=index[N]itemのサイズ
非ゼロ非対角成分総数

FEM3D101

CRS形式への変換：MAT_CON1#include <stdio.h>#include "pfem_util.h"#include "allocate.h"extern FILE* fp_log;
void MAT_CON1(){
int i,k,kk;
indexLU=(KINT*)allocate_vector(sizeof(KINT),N+1);
for(i=0;i<N+1;i++) indexLU[i]=0;
for(i=0;i<N;i++){indexLU[i+1]=indexLU[i]+INLU[i];}NPLU=indexLU[N];itemLU=(KINT*)allocate_vector(sizeof(KINT),NPLU);}

for(i=0;i<N;i++){for(k=0;k<INLU[i];k++){kk=k+indexLU[i];itemLU[kk]=IALU[i][k]-1;}}
deallocate_vector(INLU);
deallocate_vector(IALU);

itemに「0」から始まる
節点番号を記憶FEM3D

102

CRS形式への変換：MAT_CON1#include <stdio.h>#include "pfem_util.h"#include "allocate.h"extern FILE* fp_log;
void MAT_CON1(){
int i,k,kk;
indexLU=(KINT*)allocate_vector(sizeof(KINT),N+1);
for(i=0;i<N+1;i++) indexLU[i]=0;
for(i=0;i<N;i++){indexLU[i+1]=indexLU[i]+INLU[i];}NPLU=indexLU[N];itemLU=(KINT*)allocate_vector(sizeof(KINT),NPLU);}

for(i=0;i<N;i++){for(k=0;k<INLU[i];k++){kk=k+indexLU[i];itemLU[kk]=IALU[i][k]-1;}}
deallocate_vector(INLU);
deallocate_vector(IALU);

これらはもはや不要

FEM3D

103

全体処理
/**

program heat3D**/#include <stdio.h>#include <stdlib.h>
FILE* fp_log;#define GLOBAL_VALUE_DEFINE#include "pfem_util.h"//#include "solver11.h"extern void INPUT_CNTL();
extern void INPUT_GRID();
extern void MAT_CON0();
extern void MAT_CON1();
extern void MAT_ASS_MAIN();
extern void MAT_ASS_BC();
extern void SOLVE11();
extern void OUTPUT_UCD();
int main(){
INPUT_CNTL();INPUT_GRID();
MAT_CON0();
MAT_CON1();
MAT_ASS_MAIN();
MAT_ASS_BC() ;SOLVE11();}OUTPUT_UCD();FEM3D104

MAT_ASS_MAIN：全体構成
do kpn= 1, 2
do jpn= 1, 2
do ipn= 1, 2

ガウス積分点番号（ζ方向）
ガウス積分点番号（η方向）
ガウス積分点番号（ξ方向）
ガウス積分点（8個）における形状関数，
およびその「自然座標系」における微分の算出

enddo
enddo
enddo
do icel= 1, ICELTOT

要素ループ
8節点の座標から，ガウス積分点における，形状関数の「全体座標系」における微分，およびヤコビアンを算出（JACOBI）

do ie= 1, 8
do je= 1, 8

局所節点番号局所節点番号
全体節点番号：ip, jp
Aip,jpのitemLUにおけるアドレス：kk

do kpn= 1, 2
do jpn= 1, 2
do ipn= 1, 2

ガウス積分点番号（ζ方向）ie
ガウス積分点番号（η方向）
ガウス積分点番号（ξ方向）
要素積分⇒要素行列成分計算，全体行列への足しこみ

enddo
enddo
enddo
enddo
enddo
enddo

je

FEM3D105

係数行列：MAT_ASS_MAIN（1/6）#include <stdio.h>#include <math.h>#include "pfem_util.h"#include "allocate.h"extern FILE *fp_log;
extern void JACOBI();
void MAT_ASS_MAIN(){
int i,k,kk;
int ip,jp,kp;
int ipn,jpn,kpn;
int icel;
int ie,je;
int iiS,iiE;
int in1,in2,in3,in4,in5,in6,in7,in8;
double SHi;
double QP1,QM1,EP1,EM1,TP1,TM1;
double X1,X2,X3,X4,X5,X6,X7,X8;
double Y1,Y2,Y3,Y4,Y5,Y6,Y7,Y8;
double Z1,Z2,Z3,Z4,Z5,Z6,Z7,Z8;
double PNXi,PNYi,PNZi,PNXj,PNYj,PNZj;
double COND0, QV0, QVC, COEFij;
double coef;
KINT nodLOCAL[8];AMAT=(KREAL*) allocate_vector(sizeof(KREAL),NPLU);
B =(KREAL*) allocate_vector(sizeof(KREAL),N );D =(KREAL*) allocate_vector(sizeof(KREAL),N);
X =(KREAL*) allocate_vector(sizeof(KREAL),N);for(i=0;i<NPLU;i++) AMAT[i]=0.0;
for(i=0;i<N ;i++) B[i]=0.0;
for(i=0;i<N ;i++) D[i]=0.0;
for(i=0;i<N ;i++) X[i]=0.0;
WEI[0]= 1.0000000000e0;
WEI[1]= 1.0000000000e0;
POS[0]= -0.5773502692e0;
POS[1]= 0.5773502692e0;

係数行列（非零非対角成分）
右辺ベクトル
解ベクトル
係数行列（対角成分）

FEM3D106

係数行列：MAT_ASS_MAIN（1/6）#include <stdio.h>#include <math.h>#include "pfem_util.h"#include "allocate.h"extern FILE *fp_log;
extern void JACOBI();
void MAT_ASS_MAIN(){
int i,k,kk;
int ip,jp,kp;
int ipn,jpn,kpn;
int icel;
int ie,je;
int iiS,iiE;
int in1,in2,in3,in4,in5,in6,in7,in8;
double SHi;
double QP1,QM1,EP1,EM1,TP1,TM1;
double X1,X2,X3,X4,X5,X6,X7,X8;
double Y1,Y2,Y3,Y4,Y5,Y6,Y7,Y8;
double Z1,Z2,Z3,Z4,Z5,Z6,Z7,Z8;
double PNXi,PNYi,PNZi,PNXj,PNYj,PNZj;
double COND0, QV0, QVC, COEFij;
double coef;
KINT nodLOCAL[8];AMAT=(KREAL*) allocate_vector(sizeof(KREAL),NPLU);
B =(KREAL*) allocate_vector(sizeof(KREAL),N );D =(KREAL*) allocate_vector(sizeof(KREAL),N);
X =(KREAL*) allocate_vector(sizeof(KREAL),N);for(i=0;i<NPLU;i++) AMAT[i]=0.0;
for(i=0;i<N ;i++) B[i]=0.0;
for(i=0;i<N ;i++) D[i]=0.0;
for(i=0;i<N ;i++) X[i]=0.0;
WEI[0]= 1.0000000000e0;
WEI[1]= 1.0000000000e0;
POS[0]= -0.5773502692e0;
POS[1]= 0.5773502692e0;

POS： 積分点座標
WEI： 重み係数

FEM3D107

係数行列：MAT_ASS_MAIN（2/6）
/***

INIT.
PNQ
- 1st-order derivative of shape function by QSI
PNE- 1st-order derivative of shape function by ETA
PNT- 1st-order derivative of shape function by ZET***/

for(ip=0;ip<2;ip++){for(jp=0;jp<2;jp++){for(kp=0;kp<2;kp++){
QP1= 1.e0 + POS[ip];QM1= 1.e0 - POS[ip];EP1= 1.e0 + POS[jp];EM1= 1.e0 - POS[jp];TP1= 1.e0 + POS[kp];TM1= 1.e0 - POS[kp];SHAPE[ip][jp][kp][0]= O8th * QM1 * EM1 * TM1;
SHAPE[ip][jp][kp][1]= O8th * QP1 * EM1 * TM1;
SHAPE[ip][jp][kp][2]= O8th * QP1 * EP1 * TM1;
SHAPE[ip][jp][kp][3]= O8th * QM1 * EP1 * TM1;
SHAPE[ip][jp][kp][4]= O8th * QM1 * EM1 * TP1;
SHAPE[ip][jp][kp][5]= O8th * QP1 * EM1 * TP1;
SHAPE[ip][jp][kp][6]= O8th * QP1 * EP1 * TP1;
SHAPE[ip][jp][kp][7]= O8th * QM1 * EP1 * TP1;

FEM3D108

係数行列：MAT_ASS_MAIN（2/6）
/***

INIT.
PNQ
- 1st-order derivative of shape function by QSI
PNE- 1st-order derivative of shape function by ETA
PNT- 1st-order derivative of shape function by ZET***/

for(ip=0;ip<2;ip++){for(jp=0;jp<2;jp++){for(kp=0;kp<2;kp++){
QP1= 1.e0 + POS[ip];QM1= 1.e0 - POS[ip];EP1= 1.e0 + POS[jp];EM1= 1.e0 - POS[jp];TP1= 1.e0 + POS[kp];TM1= 1.e0 - POS[kp];SHAPE[ip][jp][kp][0]= O8th * QM1 * EM1 * TM1;
SHAPE[ip][jp][kp][1]= O8th * QP1 * EM1 * TM1;
SHAPE[ip][jp][kp][2]= O8th * QP1 * EP1 * TM1;
SHAPE[ip][jp][kp][3]= O8th * QM1 * EP1 * TM1;
SHAPE[ip][jp][kp][4]= O8th * QM1 * EM1 * TP1;
SHAPE[ip][jp][kp][5]= O8th * QP1 * EM1 * TP1;
SHAPE[ip][jp][kp][6]= O8th * QP1 * EP1 * TP1;
SHAPE[ip][jp][kp][7]= O8th * QM1 * EP1 * TP1;

QP1(i ) = (1 + ξ i ), QM1(i ) = (1 − ξ i )EP1( j) = (1 + η j ), EM1( j) = (1 − ηi )TP1(k ) = (1 + ζ k ), TM1(k ) = (1 − ζ k )FEM3D109

係数行列：MAT_ASS_MAIN（2/6）
/***

INIT.
PNQ
- 1st-order derivative of shape function by QSI
PNE- 1st-order derivative of shape function by ETA
PNT- 1st-order derivative of shape function by ZET***/

for(ip=0;ip<2;ip++){for(jp=0;jp<2;jp++){for(kp=0;kp<2;kp++){
QP1= 1.e0 + POS[ip];QM1= 1.e0 - POS[ip];EP1= 1.e0 + POS[jp];EM1= 1.e0 - POS[jp];TP1= 1.e0 + POS[kp];TM1= 1.e0 - POS[kp];SHAPE[ip][jp][kp][0]= O8th * QM1 * EM1 * TM1;
SHAPE[ip][jp][kp][1]= O8th * QP1 * EM1 * TM1;
SHAPE[ip][jp][kp][2]= O8th * QP1 * EP1 * TM1;
SHAPE[ip][jp][kp][3]= O8th * QM1 * EP1 * TM1;
SHAPE[ip][jp][kp][4]= O8th * QM1 * EM1 * TP1;
SHAPE[ip][jp][kp][5]= O8th * QP1 * EM1 * TP1;
SHAPE[ip][jp][kp][6]= O8th * QP1 * EP1 * TP1;
SHAPE[ip][jp][kp][7]= O8th * QM1 * EP1 * TP1;(− 1,+1,+1) 8(− 1,−1,+1) 5

7 (+ 1,+1,+1)
6 (+ 1,−1,+1)(− 1,+1,−1)
4
1(ξ ,η , ζ ) = (− 1,−1,−1)

2

3 (+ 1,+1,−1)(+ 1,−1,−1)

FEM3D110

係数行列：MAT_ASS_MAIN（2/6）
/***

INIT.
PNQ
- 1st-order derivative of shape function by QSI
PNE- 1st-order derivative of shape function by ETA
PNT- 1st-order derivative of shape function by ZET***/

for(ip=0;ip<2;ip++){for(jp=0;jp<2;jp++){for(kp=0;kp<2;kp++){
QP1= 1.e0 + POS[ip];QM1= 1.e0 - POS[ip];EP1= 1.e0 + POS[jp];EM1= 1.e0 - POS[jp];TP1= 1.e0 + POS[kp];TM1= 1.e0 - POS[kp];SHAPE[ip][jp][kp][0]= O8th * QM1 * EM1 * TM1;
SHAPE[ip][jp][kp][1]= O8th * QP1 * EM1 * TM1;
SHAPE[ip][jp][kp][2]= O8th * QP1 * EP1 * TM1;
SHAPE[ip][jp][kp][3]= O8th * QM1 * EP1 * TM1;
SHAPE[ip][jp][kp][4]= O8th * QM1 * EM1 * TP1;
SHAPE[ip][jp][kp][5]= O8th * QP1 * EM1 * TP1;
SHAPE[ip][jp][kp][6]= O8th * QP1 * EP1 * TP1;
SHAPE[ip][jp][kp][7]= O8th * QM1 * EP1 * TP1;

1
(1 − ξ )(1 − η )(1 − ζ )81
N 2 (ξ ,η , ζ ) = (1 + ξ )(1 − η )(1 − ζ )81
N 3 (ξ ,η , ζ ) = (1 + ξ )(1 + η )(1 − ζ )81
N 4 (ξ ,η , ζ ) = (1 − ξ )(1 + η )(1 − ζ )8
N1 (ξ ,η , ζ ) =1(1 − ξ )(1 − η )(1 + ζ )81
N 6 (ξ ,η , ζ ) = (1 + ξ )(1 − η )(1 + ζ )81
N 7 (ξ ,η , ζ ) = (1 + ξ )(1 + η )(1 + ζ )81
N 8 (ξ ,η , ζ ) = (1 − ξ )(1 + η )(1 + ζ )8
N 5 (ξ ,η , ζ ) =FEM3D111

係数行列：MAT_ASS_MAIN（3/6）
PNQ[jp][kp][0]= - O8th * EM1 * TM1;PNQ[jp][kp][1]= + O8th * EM1 * TM1;PNQ[jp][kp][2]= + O8th * EP1 * TM1;PNQ[jp][kp][3]= - O8th * EP1 * TM1;PNQ[jp][kp][4]= - O8th * EM1 * TP1;PNQ[jp][kp][5]= + O8th * EM1 * TP1;PNQ[jp][kp][6]= + O8th * EP1 * TP1;PNQ[jp][kp][7]= - O8th * EP1 * TP1;PNE[ip][kp][0]= - O8th * QM1 * TM1;
PNE[ip][kp][1]= - O8th * QP1 * TM1;
PNE[ip][kp][2]= + O8th * QP1 * TM1;
PNE[ip][kp][3]= + O8th * QM1 * TM1;
PNE[ip][kp][4]= - O8th * QM1 * TP1;PNE[ip][kp][5]= - O8th * QP1 * TP1;PNE[ip][kp][6]= + O8th * QP1 * TP1;PNE[ip][kp][7]= + O8th * QM1 * TP1;

}

}

}

PNT[ip][jp][0]= - O8th * QM1 * EM1;PNT[ip][jp][1]= - O8th * QP1 * EM1;PNT[ip][jp][2]= - O8th * QP1 * EP1;PNT[ip][jp][3]= - O8th * QM1 * EP1;PNT[ip][jp][4]= + O8th * QM1 * EM1;PNT[ip][jp][5]= + O8th * QP1 * EM1;PNT[ip][jp][6]= + O8th * QP1 * EP1;PNT[ip][jp][7]= + O8th * QM1 * EP1;

for( icel=0;icel< ICELTOT;icel++){COND0= COND;
in1=ICELNOD[icel][0];in2=ICELNOD[icel][1];in3=ICELNOD[icel][2];in4=ICELNOD[icel][3];in5=ICELNOD[icel][4];in6=ICELNOD[icel][5];in7=ICELNOD[icel][6];in8=ICELNOD[icel][7];∂N l(ξ = ξ i ,η = η j , ζ = ζ k )∂ξ
∂N
PNE (i, k ) = l (ξ = ξ i ,η = η j , ζ = ζ k )∂η
∂N
PNT (i, j ) = l (ξ = ξ i ,η = η j , ζ = ζ k )∂ζ
PNQ( j , k ) =∂N11(ξ i ,η j , ζ k ) = − (1 − η j )(1 − ζ k )
∂ξ
8
∂N 2
1
(ξ i ,η j , ζ k ) = + (1 − η j )(1 − ζ k )
∂ξ
8
∂N 3
1
(ξ i ,η j , ζ k ) = + (1 + η j )(1 − ζ k )
∂ξ
8
∂N 3
1
(ξ i ,η j , ζ k ) = − (1 + η j )(1 − ζ k )∂ξ
8(ξ i ,η j , ζ k ) における形状関数の一階微分

FEM3D112

係数行列：MAT_ASS_MAIN（3/6）
PNQ[jp][kp][0]= - O8th * EM1 * TM1;PNQ[jp][kp][1]= + O8th * EM1 * TM1;PNQ[jp][kp][2]= + O8th * EP1 * TM1;PNQ[jp][kp][3]= - O8th * EP1 * TM1;PNQ[jp][kp][4]= - O8th * EM1 * TP1;PNQ[jp][kp][5]= + O8th * EM1 * TP1;PNQ[jp][kp][6]= + O8th * EP1 * TP1;PNQ[jp][kp][7]= - O8th * EP1 * TP1;PNE[ip][kp][0]= - O8th * QM1 * TM1;
PNE[ip][kp][1]= - O8th * QP1 * TM1;
PNE[ip][kp][2]= + O8th * QP1 * TM1;
PNE[ip][kp][3]= + O8th * QM1 * TM1;
PNE[ip][kp][4]= - O8th * QM1 * TP1;PNE[ip][kp][5]= - O8th * QP1 * TP1;PNE[ip][kp][6]= + O8th * QP1 * TP1;PNE[ip][kp][7]= + O8th * QM1 * TP1;

}

}

}

PNT[ip][jp][0]= - O8th * QM1 * EM1;PNT[ip][jp][1]= - O8th * QP1 * EM1;PNT[ip][jp][2]= - O8th * QP1 * EP1;PNT[ip][jp][3]= - O8th * QM1 * EP1;PNT[ip][jp][4]= + O8th * QM1 * EM1;PNT[ip][jp][5]= + O8th * QP1 * EM1;PNT[ip][jp][6]= + O8th * QP1 * EP1;PNT[ip][jp][7]= + O8th * QM1 * EP1;(− 1,+1,+1) 8(− 1,−1,+1) 5

6 (+ 1,−1,+1)(− 1,+1,−1)for( icel=0;icel< ICELTOT;icel++){COND0= COND;
in1=ICELNOD[icel][0];in2=ICELNOD[icel][1];in3=ICELNOD[icel][2];in4=ICELNOD[icel][3];in5=ICELNOD[icel][4];in6=ICELNOD[icel][5];in7=ICELNOD[icel][6];in8=ICELNOD[icel][7];7 (+ 1,+1,+1)

4
1(ξ ,η , ζ ) = (− 1,−1,−1)

2

3 (+ 1,+1,−1)(+ 1,−1,−1)

FEM3D113

係数行列：MAT_ASS_MAIN（4/6）
nodLOCAL[0]= in1;
nodLOCAL[1]= in2;
nodLOCAL[2]= in3;
nodLOCAL[3]= in4;
nodLOCAL[4]= in5;
nodLOCAL[5]= in6;
nodLOCAL[6]= in7;
nodLOCAL[7]= in8;
X1=XYZ[in1-1][0];X2=XYZ[in2-1][0];X3=XYZ[in3-1][0];X4=XYZ[in4-1][0];X5=XYZ[in5-1][0];X6=XYZ[in6-1][0];X7=XYZ[in7-1][0];X8=XYZ[in8-1][0];(− 1,+1,+1) 8

8節点の節点番号(− 1,−1,+1) 5

7 (+ 1,+1,+1)
6 (+ 1,−1,+1)(− 1,+1,−1)
4
1(ξ ,η , ζ ) = (− 1,−1,−1)

Y1=XYZ[in1-1][1];Y2=XYZ[in2-1][1];Y3=XYZ[in3-1][1];Y4=XYZ[in4-1][1];Y5=XYZ[in5-1][1];Y6=XYZ[in6-1][1];Y7=XYZ[in7-1][1];Y8=XYZ[in8-1][1];QVC= O8th*(X1+X2+X3+X4+X5+X6+X7+X8+Y1+Y2+Y3+Y4+Y5+Y6+Y7+Y8);
Z1=XYZ[in1-1][2];Z2=XYZ[in2-1][2];Z3=XYZ[in3-1][2];Z4=XYZ[in4-1][2];Z5=XYZ[in5-1][2];Z6=XYZ[in6-1][2];Z7=XYZ[in7-1][2];Z8=XYZ[in8-1][2];JACOBI(DETJ, PNQ, PNE, PNT, PNX, PNY, PNZ,
X1, X2, X3, X4, X5, X6, X7, X8,
Y1, Y2, Y3, Y4, Y5, Y6, Y7, Y8, Z1, Z2, Z3, Z4, Z5, Z6, Z7, Z8);

2

3 (+ 1,+1,−1)(+ 1,−1,−1)

FEM3D114

係数行列：MAT_ASS_MAIN（4/6）
nodLOCAL[0]= in1;
nodLOCAL[1]= in2;
nodLOCAL[2]= in3;
nodLOCAL[3]= in4;
nodLOCAL[4]= in5;
nodLOCAL[5]= in6;
nodLOCAL[6]= in7;
nodLOCAL[7]= in8;
X1=XYZ[in1-1][0];X2=XYZ[in2-1][0];X3=XYZ[in3-1][0];X4=XYZ[in4-1][0];X5=XYZ[in5-1][0];X6=XYZ[in6-1][0];X7=XYZ[in7-1][0];X8=XYZ[in8-1][0];Y1=XYZ[in1-1][1];Y2=XYZ[in2-1][1];Y3=XYZ[in3-1][1];Y4=XYZ[in4-1][1];Y5=XYZ[in5-1][1];Y6=XYZ[in6-1][1];Y7=XYZ[in7-1][1];Y8=XYZ[in8-1][1];(− 1,+1,+1) 8(− 1,−1,+1) 5

6 (+ 1,−1,+1)(− 1,+1,−1)
4

8節点のX座標1(ξ ,η , ζ ) = (− 1,−1,−1)

8節点のY座標23 (+ 1,+1,−1)(+ 1,−1,−1)

座標値：
節点番号から1引く

QVC= O8th*(X1+X2+X3+X4+X5+X6+X7+X8+Y1+Y2+Y3+Y4+Y5+Y6+Y7+Y8);
Z1=XYZ[in1-1][2];Z2=XYZ[in2-1][2];Z3=XYZ[in3-1][2];Z4=XYZ[in4-1][2];Z5=XYZ[in5-1][2];Z6=XYZ[in6-1][2];Z7=XYZ[in7-1][2];Z8=XYZ[in8-1][2];

7 (+ 1,+1,+1)

8節点のZ座標

JACOBI(DETJ, PNQ, PNE, PNT, PNX, PNY, PNZ,
X1, X2, X3, X4, X5, X6, X7, X8,
Y1, Y2, Y3, Y4, Y5, Y6, Y7, Y8, Z1, Z2, Z3, Z4, Z5, Z6, Z7, Z8);

FEM3D115

係数行列：MAT_ASS_MAIN（4/6）
nodLOCAL[0]= in1;
nodLOCAL[1]= in2;
nodLOCAL[2]= in3;
nodLOCAL[3]= in4;
nodLOCAL[4]= in5;
nodLOCAL[5]= in6;
nodLOCAL[6]= in7;
nodLOCAL[7]= in8;
X1=XYZ[in1-1][0];X2=XYZ[in2-1][0];X3=XYZ[in3-1][0];X4=XYZ[in4-1][0];X5=XYZ[in5-1][0];X6=XYZ[in6-1][0];X7=XYZ[in7-1][0];X8=XYZ[in8-1][0];Y1=XYZ[in1-1][1];Y2=XYZ[in2-1][1];Y3=XYZ[in3-1][1];Y4=XYZ[in4-1][1];Y5=XYZ[in5-1][1];Y6=XYZ[in6-1][1];Y7=XYZ[in7-1][1];Y8=XYZ[in8-1][1];(− 1,+1,+1) 8

7 (+ 1,+1,+1)(− 1,−1,+1) 5

6 (+ 1,−1,+1)(− 1,+1,−1)
4

8節点のX座標1(ξ ,η , ζ ) = (− 1,−1,−1)

2

3 (+ 1,+1,−1)(+ 1,−1,−1)

座標値：
節点番号から1引く

8節点のY座標

QVC= O8th*(X1+X2+X3+X4+X5+X6+X7+X8+Y1+Y2+Y3+Y4+Y5+Y6+Y7+Y8);
Z1=XYZ[in1-1][2];Z2=XYZ[in2-1][2];Z3=XYZ[in3-1][2];Z4=XYZ[in4-1][2];Z5=XYZ[in5-1][2];Z6=XYZ[in6-1][2];Z7=XYZ[in7-1][2];Z8=XYZ[in8-1][2];∂  ∂T  ∂  ∂T  ∂  ∂T  ɺ +  λλ +  λ + Q ( x, y , z ) = 0∂x  ∂x  ∂y  ∂y  ∂z  ∂z Qɺ (x, y, z ) = QVOL xC + yC

体積当たり発熱量は位置（メッシュの中心

JACOBI(DETJ, PNQ, PNE, PNT, PNX, PNY, PNZ,
X1, X2, X3, X4, X5, X6, X7, X8,の座標 xc,yc）に依存
Y1, Y2, Y3, Y4, Y5, Y6, Y7, Y8, Z1, Z2, Z3, Z4, Z5, Z6, Z7, Z8);

FEM3D116

係数行列：MAT_ASS_MAIN（4/6）
nodLOCAL[0]= in1;
nodLOCAL[1]= in2;
nodLOCAL[2]= in3;
nodLOCAL[3]= in4;
nodLOCAL[4]= in5;
nodLOCAL[5]= in6;
nodLOCAL[6]= in7;
nodLOCAL[7]= in8;
X1=XYZ[in1-1][0];X2=XYZ[in2-1][0];X3=XYZ[in3-1][0];X4=XYZ[in4-1][0];X5=XYZ[in5-1][0];X6=XYZ[in6-1][0];X7=XYZ[in7-1][0];X8=XYZ[in8-1][0];(− 1,+1,+1) 8

7 (+ 1,+1,+1)(− 1,−1,+1) 5

6 (+ 1,−1,+1)(− 1,+1,−1)
4
1(ξ ,η , ζ ) = (− 1,−1,−1)

Y1=XYZ[in1-1][1];Y2=XYZ[in2-1][1];Y3=XYZ[in3-1][1];Y4=XYZ[in4-1][1];Y5=XYZ[in5-1][1];Y6=XYZ[in6-1][1];Y7=XYZ[in7-1][1];Y8=XYZ[in8-1][1];2

3 (+ 1,+1,−1)(+ 1,−1,−1)

座標値：
節点番号から1引く

QVC= O8th*(X1+X2+X3+X4+X5+X6+X7+X8+Y1+Y2+Y3+Y4+Y5+Y6+Y7+Y8);
Z1=XYZ[in1-1][2];Z2=XYZ[in2-1][2];Z3=XYZ[in3-1][2];Z4=XYZ[in4-1][2];Z5=XYZ[in5-1][2];Z6=XYZ[in6-1][2];Z7=XYZ[in7-1][2];Z8=XYZ[in8-1][2];∂  ∂T  ∂  ∂T  ∂  ∂T  ɺ +  λλ +  λ + Q ( x, y , z ) = 0∂x  ∂x  ∂y  ∂y  ∂z  ∂z Qɺ (x, y, z ) = QVOL xC + yC
QVC = xC + yC

JACOBI(DETJ, PNQ, PNE, PNT, PNX, PNY, PNZ,
X1, X2, X3, X4, X5, X6, X7, X8,
Y1, Y2, Y3, Y4, Y5, Y6, Y7, Y8, Z1, Z2, Z3, Z4, Z5, Z6, Z7, Z8);

FEM3D117

係数行列：MAT_ASS_MAIN（4/6）
nodLOCAL[0]= in1;
nodLOCAL[1]= in2;
nodLOCAL[2]= in3;
nodLOCAL[3]= in4;
nodLOCAL[4]= in5;
nodLOCAL[5]= in6;
nodLOCAL[6]= in7;
nodLOCAL[7]= in8;
X1=XYZ[in1-1][0];X2=XYZ[in2-1][0];X3=XYZ[in3-1][0];X4=XYZ[in4-1][0];X5=XYZ[in5-1][0];X6=XYZ[in6-1][0];X7=XYZ[in7-1][0];X8=XYZ[in8-1][0];Y1=XYZ[in1-1][1];Y2=XYZ[in2-1][1];Y3=XYZ[in3-1][1];Y4=XYZ[in4-1][1];Y5=XYZ[in5-1][1];Y6=XYZ[in6-1][1];Y7=XYZ[in7-1][1];Y8=XYZ[in8-1][1];QVC= O8th*(X1+X2+X3+X4+X5+X6+X7+X8+Y1+Y2+Y3+Y4+Y5+Y6+Y7+Y8);
Z1=XYZ[in1-1][2];Z2=XYZ[in2-1][2];Z3=XYZ[in3-1][2];Z4=XYZ[in4-1][2];Z5=XYZ[in5-1][2];Z6=XYZ[in6-1][2];Z7=XYZ[in7-1][2];Z8=XYZ[in8-1][2];JACOBI(DETJ, PNQ, PNE, PNT, PNX, PNY, PNZ,
X1, X2, X3, X4, X5, X6, X7, X8,
Y1, Y2, Y3, Y4, Y5, Y6, Y7, Y8, Z1, Z2, Z3, Z4, Z5, Z6, Z7, Z8);

FEM3D118

JACOBI（1/4）#include <stdio.h>#include <math.h>#include "precision.h"#include "allocate.h"/****** JACOBI***/
void JACOBI(KREAL DETJ[2][2][2],KREAL PNQ[2][2][8],KREAL PNE[2][2][8],KREAL PNT[2][2][8],KREAL PNX[2][2][2][8],KREAL PNY[2][2][2][8],KREAL PNZ[2][2][2][8],KREAL X1,KREAL X2,KREAL X3,KREAL X4,KREAL X5,KREAL X6,KREAL X7,KREAL X8,
KREAL Y1,KREAL Y2,KREAL Y3,KREAL Y4,KREAL Y5,KREAL Y6,KREAL Y7,KREAL Y8,
KREAL Z1,KREAL Z2,KREAL Z3,KREAL Z4,KREAL Z5,KREAL Z6,KREAL Z7,KREAL Z8){/**
calculates JACOBIAN & INVERSE JACOBIAN
dNi/dx, dNi/dy & dNi/dz**/
int ip,jp,kp;
double dXdQ,dYdQ,dZdQ,dXdE,dYdE,dZdE,dXdT,dYdT,dZdT;
double coef;
double a11,a12,a13,a21,a22,a23,a31,a32,a33;
for(ip=0;ip<2;ip++){for(jp=0;jp<2;jp++){for(kp=0;kp<2;kp++){PNX[ip][jp][kp][0]=0.0;
PNX[ip][jp][kp][1]=0.0;
PNX[ip][jp][kp][2]=0.0;
PNX[ip][jp][kp][3]=0.0;
PNX[ip][jp][kp][4]=0.0;
PNX[ip][jp][kp][5]=0.0;
PNX[ip][jp][kp][6]=0.0;
PNX[ip][jp][kp][7]=0.0;

入力 ∂N l ∂N l ∂N l  ∂ξ , ∂η , ∂ζ , ( xl , yl , zl ) (l = 1 ~ 8)出力 ∂N l ∂N l ∂N l  ∂x , ∂y , ∂z , det J

各ガウス積分点[ip][jp][kp]における値

FEM3D

119

自然座標系における偏微分（1/4）• 偏微分の公式より以下のようになる：∂N i (ξ ,η , ζ ) ∂N i ∂x ∂N i ∂y ∂N i ∂z=++
∂x ∂ξ ∂y ∂ξ ∂z ∂ξ∂ξ∂N i (ξ ,η , ζ ) ∂N i ∂x ∂N i ∂y ∂N i ∂z=++∂η∂x ∂η ∂y ∂η ∂z ∂η
∂N i (ξ ,η , ζ ) ∂N i ∂x ∂N i ∂y ∂N i ∂z
=
+
+
∂ζ
∂x ∂ζ
∂y ∂ζ
∂z ∂ζ
 ∂N i ∂N i ∂N i  ∂ξ , ∂η , ∂ζ  は定義より簡単に求められるが
 ∂N i ∂N i ∂N i  ∂x , ∂y , ∂z  を実際の計算で使用する


FEM3D

120

自然座標系における偏微分（2/4）• マトリックス表示すると： ∂N i   ∂x ∂ξ   ∂ξ

 
∂
N
 i   ∂x= ∂η   ∂η ∂N i   ∂x ∂ζ   ∂ζ

∂y
∂ξ
∂y
∂η
∂y
∂ζ

∂z   ∂N i  ∂N i  ∂x ∂ξ   ∂x  ∂N ∂z   ∂N i 
 i
[
]
=
J





∂η  ∂y ∂y 
∂z   ∂N i ∂N i  ∂z ∂ζ   ∂z  J11[J ] =  J 21
 J 31J12

J13 
J 23 
J 33 

J 22
J 32[J ]： ヤコビのマトリクス（Jacobi matrix
Jacobian）

FEM3D

121

自然座標系における偏微分（3/4）• Niの定義より簡単に求められる
∂y
∂  8
∂x
∂  8
 8 ∂N i 8 ∂N i
J11 =xi , J12 =yi ,==  N i yi  =   N i xi  = ∂ξ ∂ξ  i =1∂ξ ∂ξ  i =1
 i =1 ∂ξ i =1 ∂ξ∂z∂  8 8 ∂N i
J13 =z=
  N i zi  = ∂ξ ∂ξ  i =1
 i =1 ∂ξ
∂x
∂  8
∂y
∂  8
 8 ∂N i 8 ∂N i
J 21 =
xi , J 22 =yi ,==
  N i xi  =   N i yi  = ∂η ∂η  i =1∂η ∂η  i =1
 i =1 ∂η i =1 ∂η∂z∂  8 8 ∂N i
J 23 =zi=  N i zi  = ∂η ∂η  i =1
 i =1 ∂η
∂x
∂  8
∂y
∂  8
 8 ∂N i 8 ∂N i
J 31 =
xi , J 32 =yi ,==
  N i xi  =   N i yi  = ∂ζ ∂ζ  i =1∂ζ ∂ζ  i =1
 i =1 ∂ζ i =1 ∂ζ∂z∂  8 8 ∂N i
J 33 ==zi  N i zi  = ∂ζ ∂ζ  i =1
 i =1 ∂ζ

FEM3D122

JACOBI（2/4）

PNY[ip][jp][kp][0]=0.0;
PNY[ip][jp][kp][1]=0.0;
PNY[ip][jp][kp][2]=0.0;
PNY[ip][jp][kp][3]=0.0;
PNY[ip][jp][kp][4]=0.0;
PNY[ip][jp][kp][5]=0.0;
PNY[ip][jp][kp][6]=0.0;
PNY[ip][jp][kp][7]=0.0;
PNZ[ip][jp][kp][0]=0.0;
PNZ[ip][jp][kp][1]=0.0;
PNZ[ip][jp][kp][2]=0.0;
PNZ[ip][jp][kp][3]=0.0;
PNZ[ip][jp][kp][4]=0.0;
PNZ[ip][jp][kp][5]=0.0;
PNZ[ip][jp][kp][6]=0.0;
PNZ[ip][jp][kp][7]=0.0;
/**
**/

 J11
[J ] =  J 21
 J 31

DETERMINANT of the JACOBIAN
dXdQ = PNQ[jp][kp][0]*X1 + PNQ[jp][kp][1]*X2+ PNQ[jp][kp][2]*X3 + PNQ[jp][kp][3]*X4+ PNQ[jp][kp][4]*X5 + PNQ[jp][kp][5]*X6+ PNQ[jp][kp][6]*X7 + PNQ[jp][kp][7]*X8;
dYdQ = PNQ[jp][kp][0]*Y1 + PNQ[jp][kp][1]*Y2+ PNQ[jp][kp][2]*Y3 + PNQ[jp][kp][3]*Y4+ PNQ[jp][kp][4]*Y5 + PNQ[jp][kp][5]*Y6+ PNQ[jp][kp][6]*Y7 + PNQ[jp][kp][7]*Y8;
dZdQ = PNQ[jp][kp][0]*Z1 + PNQ[jp][kp][1]*Z2+ PNQ[jp][kp][2]*Z3 + PNQ[jp][kp][3]*Z4+ PNQ[jp][kp][4]*Z5 + PNQ[jp][kp][5]*Z6+ PNQ[jp][kp][6]*Z7 + PNQ[jp][kp][7]*Z8;
dXdE = PNE[ip][kp][0]*X1 + PNE[ip][kp][1]*X2+ PNE[ip][kp][2]*X3 + PNE[ip][kp][3]*X4+ PNE[ip][kp][4]*X5 + PNE[ip][kp][5]*X6+ PNE[ip][kp][6]*X7 + PNE[ip][kp][7]*X8;

J12
J 22
J 32

J13 
J 23 
J 33 

∂x
= J11
∂ξ
∂y
= J12
dYdQ =∂ξ∂zdZdQ == J13∂ξdXdQ =

FEM3D123

JACOBI（3/4）

/**
**/

dYdE = PNE[ip][kp][0]*Y1 + PNE[ip][kp][1]*Y2+ PNE[ip][kp][2]*Y3 + PNE[ip][kp][3]*Y4+ PNE[ip][kp][4]*Y5 + PNE[ip][kp][5]*Y6+ PNE[ip][kp][6]*Y7 + PNE[ip][kp][7]*Y8;
dZdE = PNE[ip][kp][0]*Z1 + PNE[ip][kp][1]*Z2+ PNE[ip][kp][2]*Z3 + PNE[ip][kp][3]*Z4+ PNE[ip][kp][4]*Z5 + PNE[ip][kp][5]*Z6+ PNE[ip][kp][6]*Z7 + PNE[ip][kp][7]*Z8;
dXdT = PNT[ip][jp][0]*X1 + PNT[ip][jp][1]*X2+ PNT[ip][jp][2]*X3 + PNT[ip][jp][3]*X4+ PNT[ip][jp][4]*X5 + PNT[ip][jp][5]*X6+ PNT[ip][jp][6]*X7 + PNT[ip][jp][7]*X8;
dYdT = PNT[ip][jp][0]*Y1 + PNT[ip][jp][1]*Y2+ PNT[ip][jp][2]*Y3 + PNT[ip][jp][3]*Y4+ PNT[ip][jp][4]*Y5 + PNT[ip][jp][5]*Y6+ PNT[ip][jp][6]*Y7 + PNT[ip][jp][7]*Y8;
dZdT = PNT[ip][jp][0]*Z1 + PNT[ip][jp][1]*Z2+ PNT[ip][jp][2]*Z3 + PNT[ip][jp][3]*Z4+ PNT[ip][jp][4]*Z5 + PNT[ip][jp][5]*Z6+ PNT[ip][jp][6]*Z7 + PNT[ip][jp][7]*Z8;
DETJ[ip][jp][kp]= dXdQ*(dYdE*dZdT-dZdE*dYdT) +dYdQ*(dZdE*dXdT-dXdE*dZdT) +dZdQ*(dXdE*dYdT-dYdE*dXdT);
INVERSE JACOBIAN
coef=1.0 / DETJ[ip][jp][kp];a11= coef * ( dYdE*dZdT - dZdE*dYdT );a12= coef * ( dZdQ*dYdT - dYdQ*dZdT );a13= coef * ( dYdQ*dZdE - dZdQ*dYdE );a21= coef * ( dZdE*dXdT - dXdE*dZdT );a22= coef * ( dXdQ*dZdT - dZdQ*dXdT );a23= coef * ( dZdQ*dXdE - dXdQ*dZdE );a31= coef * ( dXdE*dYdT - dYdE*dXdT );a32= coef * ( dYdQ*dXdT - dXdQ*dYdT );a33= coef * ( dXdQ*dYdE - dYdQ*dXdE );DETJ[ip][jp][kp]=fabs(DETJ[ip][jp][kp]);

 J11[J ] =  J 21
 J 31

J12
J 22
J 32

J13 
J 23 
J 33 

FEM3D

124

自然座標系における偏微分（4/4）• 従って下記のように偏微分を計算できる– ヤコビアン（3×3行列）の逆行列を求める
 ∂N i   ∂x
 ∂x   ∂ξ ∂N   ∂x
 i = ∂y   ∂η ∂N i   ∂x ∂z   ∂ζ

∂y
∂ξ
∂y
∂η
∂y
∂ζ

∂z 
∂ξ 
∂z 
∂η 
∂z 

∂ζ 

−1

 ∂N i  ∂N i  ∂ξ  ∂ξ 
 ∂N i −1  ∂N i 
[
]
=
J




η
∂
∂
η




 ∂N i  ∂N i  ∂ζ  ∂ζ 

FEM3D125

JACOBI（3/4）

/**
**/

dYdE = PNE[ip][kp][0]*Y1 + PNE[ip][kp][1]*Y2+ PNE[ip][kp][2]*Y3 + PNE[ip][kp][3]*Y4+ PNE[ip][kp][4]*Y5 + PNE[ip][kp][5]*Y6+ PNE[ip][kp][6]*Y7 + PNE[ip][kp][7]*Y8;
dZdE = PNE[ip][kp][0]*Z1 + PNE[ip][kp][1]*Z2+ PNE[ip][kp][2]*Z3 + PNE[ip][kp][3]*Z4+ PNE[ip][kp][4]*Z5 + PNE[ip][kp][5]*Z6+ PNE[ip][kp][6]*Z7 + PNE[ip][kp][7]*Z8;
dXdT = PNT[ip][jp][0]*X1 + PNT[ip][jp][1]*X2+ PNT[ip][jp][2]*X3 + PNT[ip][jp][3]*X4+ PNT[ip][jp][4]*X5 + PNT[ip][jp][5]*X6+ PNT[ip][jp][6]*X7 + PNT[ip][jp][7]*X8;
dYdT = PNT[ip][jp][0]*Y1 + PNT[ip][jp][1]*Y2+ PNT[ip][jp][2]*Y3 + PNT[ip][jp][3]*Y4+ PNT[ip][jp][4]*Y5 + PNT[ip][jp][5]*Y6+ PNT[ip][jp][6]*Y7 + PNT[ip][jp][7]*Y8;
dZdT = PNT[ip][jp][0]*Z1 + PNT[ip][jp][1]*Z2+ PNT[ip][jp][2]*Z3 + PNT[ip][jp][3]*Z4+ PNT[ip][jp][4]*Z5 + PNT[ip][jp][5]*Z6+ PNT[ip][jp][6]*Z7 + PNT[ip][jp][7]*Z8;
DETJ[ip][jp][kp]= dXdQ*(dYdE*dZdT-dZdE*dYdT) +dYdQ*(dZdE*dXdT-dXdE*dZdT) +dZdQ*(dXdE*dYdT-dYdE*dXdT);
INVERSE JACOBIAN
coef=1.0 / DETJ[ip][jp][kp];a11= coef * ( dYdE*dZdT - dZdE*dYdT );a12= coef * ( dZdQ*dYdT - dYdQ*dZdT );a13= coef * ( dYdQ*dZdE - dZdQ*dYdE );a21= coef * ( dZdE*dXdT - dXdE*dZdT );a22= coef * ( dXdQ*dZdT - dZdQ*dXdT );a23= coef * ( dZdQ*dXdE - dXdQ*dZdE );a31= coef * ( dXdE*dYdT - dYdE*dXdT );a32= coef * ( dYdQ*dXdT - dXdQ*dYdT );a33= coef * ( dXdQ*dYdE - dYdQ*dXdE );DETJ[ip][jp][kp]=fabs(DETJ[ip][jp][kp]); a11[J ]−1 = a21
 a31

a12
a22
a32

a13 
a23 
a33 

FEM3D126

JACOBI（4/4）
/**
**/

}

set the dNi/dX, dNi/dY & dNi/dZ components
PNX[ip][jp][kp][0]=a11*PNQ[jp][kp][0]+a12*PNE[ip][kp][0]+a13*PNT[ip][jp][0];PNX[ip][jp][kp][1]=a11*PNQ[jp][kp][1]+a12*PNE[ip][kp][1]+a13*PNT[ip][jp][1];PNX[ip][jp][kp][2]=a11*PNQ[jp][kp][2]+a12*PNE[ip][kp][2]+a13*PNT[ip][jp][2];PNX[ip][jp][kp][3]=a11*PNQ[jp][kp][3]+a12*PNE[ip][kp][3]+a13*PNT[ip][jp][3];PNX[ip][jp][kp][4]=a11*PNQ[jp][kp][4]+a12*PNE[ip][kp][4]+a13*PNT[ip][jp][4];PNX[ip][jp][kp][5]=a11*PNQ[jp][kp][5]+a12*PNE[ip][kp][5]+a13*PNT[ip][jp][5];PNX[ip][jp][kp][6]=a11*PNQ[jp][kp][6]+a12*PNE[ip][kp][6]+a13*PNT[ip][jp][6];PNX[ip][jp][kp][7]=a11*PNQ[jp][kp][7]+a12*PNE[ip][kp][7]+a13*PNT[ip][jp][7];PNY[ip][jp][kp][0]=a21*PNQ[jp][kp][0]+a22*PNE[ip][kp][0]+a23*PNT[ip][jp][0];PNY[ip][jp][kp][1]=a21*PNQ[jp][kp][1]+a22*PNE[ip][kp][1]+a23*PNT[ip][jp][1];PNY[ip][jp][kp][2]=a21*PNQ[jp][kp][2]+a22*PNE[ip][kp][2]+a23*PNT[ip][jp][2];PNY[ip][jp][kp][3]=a21*PNQ[jp][kp][3]+a22*PNE[ip][kp][3]+a23*PNT[ip][jp][3];PNY[ip][jp][kp][4]=a21*PNQ[jp][kp][4]+a22*PNE[ip][kp][4]+a23*PNT[ip][jp][4];PNY[ip][jp][kp][5]=a21*PNQ[jp][kp][5]+a22*PNE[ip][kp][5]+a23*PNT[ip][jp][5];PNY[ip][jp][kp][6]=a21*PNQ[jp][kp][6]+a22*PNE[ip][kp][6]+a23*PNT[ip][jp][6];PNY[ip][jp][kp][7]=a21*PNQ[jp][kp][7]+a22*PNE[ip][kp][7]+a23*PNT[ip][jp][7];PNZ[ip][jp][kp][0]=a31*PNQ[jp][kp][0]+a32*PNE[ip][kp][0]+a33*PNT[ip][jp][0];PNZ[ip][jp][kp][1]=a31*PNQ[jp][kp][1]+a32*PNE[ip][kp][1]+a33*PNT[ip][jp][1];PNZ[ip][jp][kp][2]=a31*PNQ[jp][kp][2]+a32*PNE[ip][kp][2]+a33*PNT[ip][jp][2];PNZ[ip][jp][kp][3]=a31*PNQ[jp][kp][3]+a32*PNE[ip][kp][3]+a33*PNT[ip][jp][3];PNZ[ip][jp][kp][4]=a31*PNQ[jp][kp][4]+a32*PNE[ip][kp][4]+a33*PNT[ip][jp][4];PNZ[ip][jp][kp][5]=a31*PNQ[jp][kp][5]+a32*PNE[ip][kp][5]+a33*PNT[ip][jp][5];PNZ[ip][jp][kp][6]=a31*PNQ[jp][kp][6]+a32*PNE[ip][kp][6]+a33*PNT[ip][jp][6];PNZ[ip][jp][kp][7]=a31*PNQ[jp][kp][7]+a32*PNE[ip][kp][7]+a33*PNT[ip][jp][7];}−1} ∂N i 


}
x
y
z
∂
∂
∂
 ∂N  ∂x   ∂ξ ∂N   ∂x
 i = ∂y   ∂η ∂N i   ∂x ∂z   ∂ζ
i

∂ξ
∂y
∂η
∂y
∂ζ

∂ξ 
∂z 
∂η 
∂z 

∂ζ 

 ∂ξ   a11 ∂N i   = a21 ∂η   a ∂N i   31 ∂ζ 

a12
a22
a32

 ∂N i  ∂ξ 
a13   ∂N i 
a23  ∂ηa33   ∂N  i ∂ζ 

FEM3D127

係数行列：MAT_ASS_MAIN（5/6）/**

CONSTRUCT the GLOBAL MATRIX**/

for(ie=0;ie<8;ie++){ip=nodLOCAL[ie];全体行列の非対角成分

for(je=0;je<8;je++){jp=nodLOCAL[je];Aip , jpkk=-1;
if( jp != ip ){iiS=indexLU[ip-1];
iiE=indexLU[ip ];for( k=iiS;k<iiE;k++){if( itemLU[k] == jp-1 ){kk=k;break;}}}(− 1,+1,+1) 8(− 1,−1,+1) 5

6 (+ 1,−1,+1)
4(ξ ,η , ζ ) = (− 1,−1,−1)

ip= nodLOCAL[ie]jp= nodLOCAL[je]7 (+ 1,+1,+1)(− 1,+1,−1)1kk：itemにおけるアドレス2

3 (+ 1,+1,−1)(+ 1,−1,−1)

1から始まる節点番号

FEM3D128

要素マトリクス：8×8行列
j[k ] (i, j = 1…8)ij

i(− 1,+1,+1) 8(− 1,−1,+1) 5

7 (+ 1,+1,+1)
6 (+ 1,−1,+1)(− 1,+1,−1)
4
1(ξ ,η , ζ ) = (− 1,−1,−1)

2

3 (+ 1,+1,−1)(+ 1,−1,−1)

FEM3D129

係数行列：MAT_ASS_MAIN（5/6）/**

CONSTRUCT the GLOBAL MATRIX**/

for(ie=0;ie<8;ie++){ip=nodLOCAL[ie];for(je=0;je<8;je++){jp=nodLOCAL[je];要素マトリクス（ie~je）
全体マトリクス（ip~jp）の関係

kk=-1;
if( jp != ip ){iiS=indexLU[ip-1];
iiE=indexLU[ip ];for( k=iiS;k<iiE;k++){if( itemLU[k] == jp-1 ){kk=k;break;}}}(− 1,+1,+1) 8(− 1,−1,+1) 5

7 (+ 1,+1,+1)
6 (+ 1,−1,+1)(− 1,+1,−1)
4
1(ξ ,η , ζ ) = (− 1,−1,−1)kk：itemLUにおけるアドレス2

3 (+ 1,+1,−1)(+ 1,−1,−1)

FEM3D130

係数行列：MAT_ASS_MAIN（6/6）
QV0= 0.e0;
COEFij= 0.e0;for(kpn=0;kpn<2;kpn++){for(jpn=0;jpn<2;jpn++){for(ipn=0;ipn<2;ipn++){coef= fabs(DETJ[ipn][jpn][kpn])*WEI[ipn]*WEI[jpn]*WEI[kpn];
PNXi= PNX[ipn][jpn][kpn][ie];PNYi= PNY[ipn][jpn][kpn][ie];PNZi= PNZ[ipn][jpn][kpn][ie];PNXj= PNX[ipn][jpn][kpn][je];PNYj= PNY[ipn][jpn][kpn][je];PNZj= PNZ[ipn][jpn][kpn][je];COEFij+= coef*COND0*(PNXi*PNXj+PNYi*PNYj+PNZi*PNZj);
SHi= SHAPE[ipn][jpn][kpn][ie];QV0+= SHi * QVOL * coef;

}

}

}

if (jp==ip) {D[ip-1]+= COEFij;
B[ip-1]+= QV0*QVC;}
if (jp != ip) {AMAT[kk]+= COEFij;
}

}

}

}

}

 ∂N i ∂N j
∂N i ∂N j∂N i ∂N j 
−1 −1 −1 λ ∂x ∂x + λ ∂y ∂y + λ ∂z ∂z  det J dξdηdζ

+1 +1 +1

FEM3D131

係数行列：MAT_ASS_MAIN（6/6）
QV0= 0.e0;
COEFij= 0.e0;for(kpn=0;kpn<2;kpn++){for(jpn=0;jpn<2;jpn++){for(ipn=0;ipn<2;ipn++){coef= fabs(DETJ[ipn][jpn][kpn])*WEI[ipn]*WEI[jpn]*WEI[kpn];
PNXi= PNX[ipn][jpn][kpn][ie];PNYi= PNY[ipn][jpn][kpn][ie];PNZi= PNZ[ipn][jpn][kpn][ie];PNXj= PNX[ipn][jpn][kpn][je];PNYj= PNY[ipn][jpn][kpn][je];PNZj= PNZ[ipn][jpn][kpn][je];COEFij+= coef*COND0*(PNXi*PNXj+PNYi*PNYj+PNZi*PNZj);
SHi= SHAPE[ipn][jpn][kpn][ie];QV0+= SHi * QVOL * coef;

}

}

}

}

}

}

  f (ξ ,η , ζ ) dξdηdζ

I=

−1 −1 −1

if (jp==ip) {D[ip-1]+= COEFij;
B[ip-1]+= QV0*QVC;}
if (jp != ip) {AMAT[kk]+= COEFij;}}+1 +1 +1

L

M

N

i =1

j =1

k =1

[

=    Wi ⋅ W j ⋅Wk ⋅ f (ξ i ,η j , ζ k ) ∂N i ∂N j
∂N i ∂N j∂N i ∂N j 
−1 −1 −1 λ ∂x ∂x + λ ∂y ∂y + λ ∂z ∂z  det J dξdηdζ

+1 +1 +1]

FEM3D132

MAT_ASS_MAIN (6/6)
for(kpn=0;kpn<2;kpn++){for(jpn=0;jpn<2;jpn++){for(ipn=0;ipn<2;ipn++){coef= fabs(DETJ[ipn][jpn][kpn])*WEI[ipn]*WEI[jpn]*WEI[kpn];
PNXi= PNX[ipn][jpn][kpn][ie]; coef = Wi ⋅ W j ⋅ Wk ⋅ det J (ξ i ,η j , ζ k )PNYi= PNY[ipn][jpn][kpn][ie];+1 +1 +1
PNZi= PNZ[ipn][jpn][kpn][ie];I =    f (ξ ,η , ζ ) dξdηdζ
PNXj= PNX[ipn][jpn][kpn][je];−1 −1 −1
PNYj= PNY[ipn][jpn][kpn][je];LMN
PNZj= PNZ[ipn][jpn][kpn][je]; =    Wi ⋅ W j ⋅ Wk ⋅ f (ξ i ,η j , ζ k )
i =1

j =1

k =1

[

COEFij+= coef*COND0*(PNXi*PNXj+PNYi*PNYj+PNZi*PNZj);

}

}

}

 ∂N i ∂N j
∂N i ∂N j∂N i ∂N j 
+
+
λ
λ
λ
 det J dξdηdζ
−1 −1 −1  ∂x ∂x∂y ∂y∂z ∂z +1 +1 +1]

FEM3D133

MAT_ASS_MAIN (6/6)
for(kpn=0;kpn<2;kpn++){for(jpn=0;jpn<2;jpn++){for(ipn=0;ipn<2;ipn++){coef= fabs(DETJ[ipn][jpn][kpn])*WEI[ipn]*WEI[jpn]*WEI[kpn];
PNXi= PNX[ipn][jpn][kpn][ie]; coef = Wi ⋅ W j ⋅ Wk ⋅ det J (ξ i ,η j , ζ k )PNYi= PNY[ipn][jpn][kpn][ie];+1 +1 +1
PNZi= PNZ[ipn][jpn][kpn][ie];I =    f (ξ ,η , ζ ) dξdηdζ
PNXj= PNX[ipn][jpn][kpn][je];−1 −1 −1
PNYj= PNY[ipn][jpn][kpn][je];LMN
PNZj= PNZ[ipn][jpn][kpn][je]; =    Wi ⋅ W j ⋅ Wk ⋅ f (ξ i ,η j , ζ k )
i =1

j =1

k =1

[

COEFij+= coef*COND0*(PNXi*PNXj+PNYi*PNYj+PNZi*PNZj);

}

}

}

 ∂N i ∂N j
∂N i ∂N j∂N i ∂N j 
+
+
λ
λ
λ
 det J dξdηdζ
−1 −1 −1  ∂x ∂x∂y ∂y∂z ∂z +1 +1 +1]

FEM3D134

係数行列：MAT_ASS_MAIN（6/6）
QV0= 0.e0;
COEFij= 0.e0;for(kpn=0;kpn<2;kpn++){for(jpn=0;jpn<2;jpn++){for(ipn=0;ipn<2;ipn++){coef= fabs(DETJ[ipn][jpn][kpn])*WEI[ipn]*WEI[jpn]*WEI[kpn];
PNXi= PNX[ipn][jpn][kpn][ie];PNYi= PNY[ipn][jpn][kpn][ie];PNZi= PNZ[ipn][jpn][kpn][ie];PNXj= PNX[ipn][jpn][kpn][je];PNYj= PNY[ipn][jpn][kpn][je];PNZj= PNZ[ipn][jpn][kpn][je];COEFij+= coef*COND0*(PNXi*PNXj+PNYi*PNYj+PNZi*PNZj);}}

if (jp==ip) {D[ip-1]+= COEFij;
B[ip-1]+= QV0*QVC;}
if (jp != ip) {AMAT[kk]+= COEFij;
}

}

}

}

}

ij

j

SHi= SHAPE[ipn][jpn][kpn][ie];QV0+= SHi * QVOL * coef;}[k ] (i, j = 1…8)i

FEM3D135

係数行列：MAT_ASS_MAIN（6/6）
QV0= 0.e0;
COEFij= 0.e0;for(kpn=0;kpn<2;kpn++){for(jpn=0;jpn<2;jpn++){for(ipn=0;ipn<2;ipn++){coef= fabs(DETJ[ipn][jpn][kpn])*WEI[ipn]*WEI[jpn]*WEI[kpn];
PNXi= PNX[ipn][jpn][kpn][ie];PNYi= PNY[ipn][jpn][kpn][ie];PNZi= PNZ[ipn][jpn][kpn][ie];PNXj= PNX[ipn][jpn][kpn][je];PNYj= PNY[ipn][jpn][kpn][je];PNZj= PNZ[ipn][jpn][kpn][je];[k ](e) {φ }(e) = { f }(e)[ f ](e) =  Qɺ [N ]T dV

COEFij+= coef*COND0*(PNXi*PNXj+PNYi*PNYj+PNZi*PNZj);
SHi= SHAPE[ipn][jpn][kpn][ie];QV0+= SHi * QVOL * coef;

}

}

}

if (jp==ip) {D[ip-1]+= COEFij;
B[ip-1]+= QV0*QVC;}
if (jp != ip) {AMAT[kk]+= COEFij;
}

}

}

}

}

V

Qɺ ( x, y, z ) = QVOL xC + yC

QVC = xC + yC

QV 0 =  QVOL[N ] dVTV

[ f ](e) = QV 0 ⋅ QVC

FEM3D136

MAT_ASS_BC：全体構成
do i= 1, N

節点ループ
（ディリクレ）境界条件を設定する節点をマーク（IWKX）

enddo
do i= 1, N 節点ループ
if (IWKX(i,1).eq.1) then

マークされた節点だったら対応する右辺ベクトル（B）の成分，対角成分（D）の成分の修正（行・列）Z

do k= index(i-1)+1, index(i)対応する非零非対角成分（AMAT）の成分の修正（行）

enddo
endif
enddo

T=0@Z=zmax

do i= 1, N 節点ループ
do k= indexLU (i-1)+1, index (i)if (IWKX(item (k),1).eq.1) then

NZ

対応する非零非対角成分の節点がマークされていたら対応する右辺ベクトル，非零非対角成分（AMAT）の成分の修正（列）

endif
enddo
enddo

Y
NX
NY
X

FEM3D137

境界条件：MAT_ASS_BC（1/2）#include <stdio.h>#include <stdlib.h>#include <string.h>#include "pfem_util.h"#include "allocate.h"extern FILE *fp_log;
void MAT_ASS_BC(){
int i,j,k,in,ib,ib0,icel;
int in1,in2,in3,in4,in5,in6,in7,in8;
int iq1,iq2,iq3,iq4,iq5,iq6,iq7,iq8;
int iS,iE;
double STRESS,VAL;

節点グループ名が「Zmax」である
節点in（1から始まる）において：

IWKX=(KINT**) allocate_matrix(sizeof(KINT),N,2);
for(i=0;i<N;i++) for(j=0;j<2;j++) IWKX[i][j]=0;/****/

Z=ZmaxIWKX[in-1][0]= 1とする

for(in=0;in<N;in++) IWKX[in][0]=0;
ib0=-1;for( ib0=0;ib0<NODGRPtot;ib0++){if( strcmp(NODGRP_NAME[ib0].name,"Zmax") == 0 ) break;}for( ib=NODGRP_INDEX[ib0];ib<NODGRP_INDEX[ib0+1];ib++){in=NODGRP_ITEM[ib];IWKX[in-1][0]=1;}

FEM3D138

境界条件：MAT_ASS_BC（2/2）
for(in=0;in<N;in++){if( IWKX[in][0] == 1 ){B[in]= 0.e0;
D[in]= 1.e0;
for(k=indexLU[in];k<indexLU[in+1];k++){AMAT[k]= 0.e0;
}
}
}

}

for(in=0;in<N;in++){for(k=indexLU[in];k<indexLU[in+1];k++){if (IWKX[itemLU[k]][0] == 1 ) {AMAT[k]= 0.e0;}}}

FEM3D139

対象とする問題：一次元熱伝導問題
体積当たり一様発熱 Qɺ∂  ∂T  ɺ
+Q = 0λ∂x  ∂x 
x= xmaxx=0 (xmin)• 一様な：断面積A，熱伝導率λ• 体積当たり一様発熱（時間当たり）〔QL-3T-1〕 Qɺ
• 境界条件– x=0 ：Τ=0–

– x=xmax：（固定）∂T= 0（断熱）∂x

復習：一次元問題

FEM1D140

x=0で成立する方程式
T0=0
体積当たり一様発熱 Qɺ∂  ∂T  ɺ
λ+Q = 0∂x  ∂x 
x= xmaxx=0 (xmin)• 一様な：断面積A，熱伝導率λ• 体積当たり一様発熱（時間当たり）〔QL-3T-1〕 Qɺ
• 境界条件– x=0 ：Τ=0–

– x=xmax：（固定）∂T= 0（断熱）∂x

復習：一次元問題

FEM1D

141

プログラム：1d.c（6/6）
第一種境界条件@x=0/*// +---------------------+// | BOUNDARY conditions |// +---------------------+*//* X=Xmin */i=0;
jS= Index[i];AMat[jS]= 0.0;
Diag[i ]= 1.0;
Rhs [i ]= 0.0;T0=0
対角成分=1，右辺=0，非対角成分=0for(k=0;k<NPLU;k++){if(Item[k]==0){AMat[k]=0.0;}}

復習：一次元問題

FEM1D

142

プログラム：1d.c（6/6）
第一種境界条件@x=0/*// +---------------------+// | BOUNDARY conditions |// +---------------------+*//* X=Xmin */i=0;
jS= Index[i];AMat[jS]= 0.0;
Diag[i ]= 1.0;
Rhs [i ]= 0.0;T0=0
対角成分=1，右辺=0，非対角成分=0

ゼロクリア

for(k=0;k<NPLU;k++){if(Item[k]==0){AMat[k]=0.0;}}

復習：一次元問題

FEM1D

143

プログラム：1d.c（6/6）
第一種境界条件@x=0/*// +---------------------+// | BOUNDARY conditions |// +---------------------+*//* X=Xmin */i=0;
jS= Index[i];AMat[jS]= 0.0;
Diag[i ]= 1.0;
Rhs [i ]= 0.0;T0=0
対角成分=1，右辺=0，非対角成分=0
消去，ゼロクリアfor(k=0;k<NPLU;k++){if(Item[k]==0){AMat[k]=0.0;}}行列の対称性を保つため，第一種境界条
件を適用している節点に対応する「列」を，右辺に移項して消去する（今の場合は非対角成分を0にするだけで良い）

復習：一次元問題

FEM1D144

第一種境界条件がT≠0の場合/*// +---------------------+// | BOUNDARY conditions |// +---------------------+*//* X=Xmin */i=0;
jS= Index[i];AMat[jS]= 0.0;
Diag[i ]= 1.0;
Rhs [i ]= PHImin;行列の対称性を保つため，第一種境界条
件を適用している節点に対応する「列」を，右辺に移項して消去する

Diag jφ j +Index[ j +1]−1

 Amat φk = Index[ j ]kItem[ k ]= Rhs jfor(j=1;i<N;i++){for (k=Index[j];k<Index[j+1];k++){if(Item[k]==0){
Rhs [j]= Rhs[j] – Amat[k]*PHImin;
AMat[k]= 0.0;}}

復習：一次元問題

FEM1D145

第一種境界条件がT≠0の場合/*// +---------------------+// | BOUNDARY conditions |// +---------------------+*//* X=Xmin */i=0;
jS= Index[i];AMat[jS]= 0.0;
Diag[i ]= 1.0;
Rhs [i ]= PHImin;行列の対称性を保つため，第一種境界条
件を適用している節点に対応する「列」を，右辺に移項して消去する

Diag jφ j +Index[ j +1]−1

 Amat φkk = Index[ j ], k ≠ k s

Item[ k ]= Rhs j − Amatk s φItem[ k s ]= Rhs j − Amatk s φminwhere Item[k s ] = 0

for(j=1;i<N;i++){for (k=Index[j];k<Index[j+1];k++){if(Item[k]==0){
Rhs [j]= Rhs[j] – Amat[k]*PHImin;
AMat[k]= 0.0;}}

復習：一次元問題

FEM3D146

境界条件：MAT_ASS_BC（2/2）
for(in=0;in<N;in++){IWKX[in-1][0]=1となる節点に対してif( IWKX[in][0] == 1 ){B[in]= 0.e0;
対角成分=1，右辺=0，非零対角成分=0D[in]= 1.e0;
for(k=indexLU[in];k<indexLU[in+1];k++){AMAT[k]= 0.e0;
ゼロクリア
}
}
}

}

for(in=0;in<N;in++){for(k=indexLU[in];k<indexLU[in+1];k++){if (IWKX[itemLU[k]][0] == 1 ) {AMAT[k]= 0.e0;}}}

ここでやっていることも一次元の時と全く変わらないFEM3D147

境界条件：MAT_ASS_BC（2/2）
for(in=0;in<N;in++){if( IWKX[in][0] == 1 ){B[in]= 0.e0;
D[in]= 1.e0;
for(k=indexLU[in];k<indexLU[in+1];k++){AMAT[k]= 0.e0;
IWKX[in-1][0]=1となる節点を非零非対}角成分として有している節点に対して，}}
右辺へ移項，当該非零非対角成分=0for(in=0;in<N;in++){for(k=indexLU[in];k<indexLU[in+1];k++){if (IWKX[itemLU[k]][0] == 1 ) {AMAT[k]= 0.e0;
}
}
}
}
消去，
ゼロクリア

ここでやっていることも一次元の時と全く変わらないFEM3D148

test1

input_cntl

メインプログラム制御データ入力

input_grid

find_node

メッシュファイル入力

節点探索

mat_con0mSORT

行列コネクティビティ生成

ソート

mat_con1
行列コネクティビティ生成

mat_ass_main

jacobi係数行列生成ヤコビアン計算

mat_ass_bc
境界条件処理三次元熱伝導解析コードheat3Dの構成

solve11cg

線形ソルバー制御

CG法計算

output_ucd
可視化処理FEM3D

149

全体処理
/**

program heat3D**/#include <stdio.h>#include <stdlib.h>
FILE* fp_log;#define GLOBAL_VALUE_DEFINE#include "pfem_util.h"//#include "solver11.h"extern void INPUT_CNTL();
extern void INPUT_GRID();
extern void MAT_CON0();
extern void MAT_CON1();
extern void MAT_ASS_MAIN();
extern void MAT_ASS_BC();
extern void SOLVE11();
extern void OUTPUT_UCD();
int main(){
INPUT_CNTL();INPUT_GRID();
MAT_CON0();
MAT_CON1();
MAT_ASS_MAIN();
MAT_ASS_BC() ;SOLVE11();}OUTPUT_UCD();FEM3D150

SOLVE11#include <stdio.h>#include <string.h>#include <math.h>#include "pfem_util.h"#include "allocate.h"extern FILE *fp_log;
extern void CG();
void SOLVE11(){
int i,j,k,ii,L;int ERROR, ICFLAG=0;
CHAR_LENGTH BUF;/**+------------+| PARAMETERs |+------------+

**/
ITER
RESID

= pfemIarray[0];= pfemRarray[0];CG法の最大反復回数
CG法の収束判定値/**+------------------+| ITERATIVE solver |+------------------+**/CG (N,NPLU, D, AMAT, indexLU, itemLU, B, X, RESID, ITER, &ERROR);
ITERactual= ITER;}FEM3D151

前処理付き共役勾配法
Preconditioned Conjugate Gradient Method （CG）
Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)ρi-1= r(i-1) z(i-1)
if i=1p(1)= z(0)else
βi-1= ρi-1/ρi-2p(i)= z(i-1) + βi-1 p(i-1)
endifq(i)= [A]p(i)αi = ρi-1/p(i)q(i)x(i)= x(i-1) + αip(i)r(i)= r(i-1) - αiq(i)check convergence |r|end

前処理：対角スケーリング

FEM3D152

対角スケーリング，点ヤコビ前処理• 前処理行列として，もとの行列の対角成分のみを取り出した行列を前処理行列 [M] とする。
– 対角スケーリング，点ヤコビ（point-Jacobi）前処理

 D1
0

[M ] =  ...

0
 0

0

...

D2

0
0

...
0
0

DN −1
...
0

0 
0 
... 

0 
DN 

• solve [M]z(i-1)= r(i-1)という場合に逆行列を簡単に求めることができる。

FEM3D

153

CGソルバー（1/6）#include <stdio.h>#include <math.h>#include "precision.h"#include "allocate.h"extern FILE *fp_log;
void

CG

(
KINT N,KINT NPLU,KREAL D[],KREAL AMAT[],KINT indexLU[], KINT itemLU[],KREAL B[],KREAL X[],KREAL RESID,KINT ITER, KINT *ERROR){
int i,j,k;
int ieL,isL,ieU,isU;
double WVAL;
double BNRM20,BNRM2,DNRM20,DNRM2;
double S1_TIME,E1_TIME;
double ALPHA,BETA;
double C1,C10,RHO,RHO0,RHO1;
int
iterPRE;KREAL **WW;KINT R=0,Z=1,Q=1,P=2,DD=3;
KINT MAXIT;KREAL TOL;FEM3D154

CG法向け変数表
変数名

種別

サイズ

内 容

I/O

N,NP

I

I

節点数

NPLU

I

O

非零非対角成分数

D

R

[N]

O

全体マトリクス：対角成分

B,X

R

[N]

O

右辺ベクトルAMATR[NPLU]

O

全体マトリクス：非零非対角成分

indexLU

I

[N+1]

O

全体マトリクス：各行の非零非対角成分数

itemLUI[NPLU]

O

全体マトリクス：列番号

ITER

I

I/O

CG法反復回数（I：上限，O：実際の反復回
数）

RESID

R

I/O

CG法反復打ち切り残差（I；基準値，O：最終残差ノルム）

MAXIT

I

-

CG法最大反復回数

TOL

R

-

CG法反復打ち切り残差基準値

WW

R

-

CG法ワーク配列

P,Q,R,Z,DD

I-

WWの添字（ベクトルID）[4][N]

FEM3D

155

CGソルバー（1/6）#include <stdio.h>#include <math.h>#include "precision.h"#include "allocate.h"extern FILE *fp_log;(0)= b-[A]x(0)Computer
WW[i][0]= WW[i][R] ⇒ {r}for i= 1, 2, …
WW[i][1]= WW[i][Z] ⇒ {z}solve [M]z(i-1)= r(i-1)ρi-1= r(i-1) z(i-1)
WW[i][1]=WW[i][Q] ⇒ {q}
void
CG (
if i=1KINT N,KINT NPLU,KREAL D[],WW[i][2]=WW[i][P]⇒{p}
KREAL AMAT[],KINT indexLU[], KINT itemLU[],
p(1)= z(0)KREAL B[],KREAL
X[],KREAL RESID,KINT
WW[i][3]=WW[i][DD]⇒
1/{D}ITER, KINT *ERROR)
else{
βi-1= ρi-1/ρi-2
int i,j,k;
int ieL,isL,ieU,isU;
p(i)= z(i-1) + βi-1 p(i-1)
double WVAL;
endif
double BNRM20,BNRM2,DNRM20,DNRM2;
q(i)= [A]p(i)double S1_TIME,E1_TIME;
αi = ρi-1/p(i)q(i)double ALPHA,BETA;
double C1,C10,RHO,RHO0,RHO1;x(i)= x(i-1) + αip(i)intiterPRE;r(i)= r(i-1) - αiq(i)check convergence |r|
KREAL **WW;
end

KINT R=0,Z=1,Q=1,P=2,DD=3;
KINT MAXIT;KREAL TOL;FEM3D

156

CGソルバー（2/6）
WW=(KREAL**) allocate_matrix(sizeof(KREAL),4,N);
MAXIT = ITER;TOL= RESID;for(i=0;i<N;i++){X[i]=0.0;}
for(i=0;i<N;i++) for(j=0;j<4;j++) WW[j][i]=0.0;/**+-----------------------+| {r0}= {b} - [A]{xini} |+-----------------------+**/
for(j=0;j<N;j++){WW[DD][j]= 1.0/D[j];WVAL= B[j] - D[j]*X[j];for( k=indexLU[j];k<indexLU[j+1];k++){i= itemLU[k];WVAL+= -AMAT[k]*X[i];}WW[R][j]= WVAL;}

WW[i][0]= WW[i][R] ⇒ {r}WW[i][1]= WW[i][Z] ⇒ {z}WW[i][1]= WW[i][Q] ⇒ {q}WW[i][2]= WW[i][P] ⇒ {p}WW[i][3]= WW[i][DD]⇒ 1/{D}対角成分の逆数（前処理用）
その都度，除算をすると効率が悪いため，予め配列に格納

FEM3D

157

CGソルバー（2/6）
WW=(KREAL**) allocate_matrix(sizeof(KREAL),4,N);

Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)
for(i=0;i<N;i++){(i-1) z(i-1)
ρ
=
r
i-1
X[i]=0.0;
if i=1}for(i=0;i<N;i++) for(j=0;j<4;j++) WW[j][i]=0.0;p(1)= z(0)/**else+-----------------------+
βi-1= ρi-1/ρi-2| {r0}= {b} - [A]{xini} |p(i)= z(i-1) + βi-1 p(i-1)+-----------------------+
**/
endif
for(j=0;j<N;j++){q(i)= [A]p(i)WW[DD][j]= 1.0/D[j];αi = ρi-1/p(i)q(i)WVAL= B[j] - D[j]*X[j];x(i)= x(i-1) + αip(i)for( k=indexLU[j];k<indexLU[j+1];k++){r(i)= r(i-1) - αiq(i)i= itemLU[k];check convergence |r|
WVAL+= -AMAT[k]*X[i];end}
MAXIT = ITER;TOL= RESID;

WW[R][j]= WVAL;}FEM3D

158

CGソルバー（3/6）
BNRM20= 0.e0;
for(i=0;i<N;i++){
BNRM20+= B[i]*B[i];}BNRM2= BNRM20;

BNRM2=|b|2
あとで収束判定に使用

if (BNRM2 == 0.e0) BNRM2= 1.e0;
ITER = 0;
for( ITER=1;ITER<= MAXIT;ITER++){/*************************************************** Conjugate Gradient Iteration**//**+----------------+| {z}= [Minv]{r} |+----------------+**/
for(i=0;i<N;i++){WW[Z][i]= WW[DD][i]*WW[R][i];}FEM3D

159

CGソルバー（3/6）
BNRM20= 0.e0;
for(i=0;i<N;i++){
BNRM20+= B[i]*B[i];}Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)
BNRM2= BNRM20;
ρi-1= r(i-1) z(i-1)
if i=1
if (BNRM2 == 0.e0) BNRM2= 1.e0;p(1)= z(0)elseITER = 0;
βi-1= ρi-1/ρi-2for( ITER=1;ITER<= MAXIT;ITER++){p(i)= z(i-1) + βi-1 p(i-1)/**
endifGradient Iteration************************************************* Conjugate**/q(i)= [A]p(i)αi = ρi-1/p(i)q(i)/**x(i)= x(i-1) + αip(i)+----------------+
r(i)= r(i-1) - αiq(i)| {z}= [Minv]{r} |+----------------+
check convergence |r|**/end
for(i=0;i<N;i++){WW[Z][i]= WW[DD][i]*WW[R][i];}FEM3D160

CGソルバー（4/6）/**+---------------+| {RHO}= {r}{z} |+---------------+**/
RHO0= 0.e0;
for(i=0;i<N;i++){
RHO0+= WW[R][i]*WW[Z][i];}RHO= RHO0;/**+-----------------------------+| {p} = {z} if
ITER=1|| BETA= RHO / RHO1 otherwise |+-----------------------------+**/
if( ITER == 1 ){for(i=0;i<N;i++){
WW[P][i]=WW[Z][i];}}else{
BETA= RHO / RHO1;
for(i=0;i<N;i++){
WW[P][i]=WW[Z][i] + BETA*WW[P][i];}}Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)ρi-1= r(i-1) z(i-1)
if i=1p(1)= z(0)else
βi-1= ρi-1/ρi-2p(i)= z(i-1) + βi-1 p(i-1)
endifq(i)= [A]p(i)αi = ρi-1/p(i)q(i)x(i)= x(i-1) + αip(i)r(i)= r(i-1) - αiq(i)check convergence |r|end

FEM3D

161

CGソルバー（5/6）/**+-------------+| {q}= [A]{p} |+-------------+**/
for( j=0;j<N;j++){
WVAL= D[j] * WW[P][j];for(k=indexLU[j];k<indexLU[j+1];k++){i=itemLU[k];WVAL+= AMAT[k] * WW[P][i];}WW[Q][j]=WVAL;}/**+---------------------+| ALPHA= RHO / {p}{q} |+---------------------+**/
C10= 0.e0;
for(i=0;i<N;i++){C10+=WW[P][i]*WW[Q][i];}C1=C10;
ALPHA= RHO / C1;

Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)ρi-1= r(i-1) z(i-1)
if i=1p(1)= z(0)else
βi-1= ρi-1/ρi-2p(i)= z(i-1) + βi-1 p(i-1)
endifq(i)= [A]p(i)αi = ρi-1/p(i)q(i)x(i)= x(i-1) + αip(i)r(i)= r(i-1) - αiq(i)check convergence |r|end

FEM3D

162

CGソルバー（6/6）/**+----------------------+| {x}= {x} + ALPHA*{p} || {r}= {r} - ALPHA*{q} |+----------------------+**/
for(i=0;i<N;i++){X [i]+= ALPHA *WW[P][i];WW[R][i]+= -ALPHA *WW[Q][i];}DNRM20= 0.e0;
for(i=0;i<N;i++){DNRM20+=WW[R][i]*WW[R][i];}DNRM2= DNRM20;
RESID= sqrt(DNRM2/BNRM2);
if ( RESID <= TOL) break;if ( ITER == MAXIT ) *ERROR= -300;
RHO1 = RHO ;}Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)ρi-1= r(i-1) z(i-1)
if i=1p(1)= z(0)else
βi-1= ρi-1/ρi-2p(i)= z(i-1) + βi-1 p(i-1)
endifq(i)= [A]p(i)αi = ρi-1/p(i)q(i)x(i)= x(i-1) + αip(i)r(i)= r(i-1) - αiq(i)check convergence |r|end

FEM3D

163

CGソルバー（6/6）/**+----------------------+| {x}= {x} + ALPHA*{p} || {r}= {r} - ALPHA*{q} |+----------------------+**/
for(i=0;i<N;i++){X [i]+= ALPHA *WW[P][i];WW[R][i]+= -ALPHA *WW[Q][i];}DNRM20= 0.e0;
for(i=0;i<N;i++){DNRM20+=WW[R][i]*WW[R][i];}DNRM2= DNRM20;
RESID= sqrt(DNRM2/BNRM2);
if ( RESID <= TOL) break;if ( ITER == MAXIT ) *ERROR= -300;
RHO1 = RHO ;}Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)ρi-1= r(i-1) z(i-1)
if i=1p(1)= z(0)else
βi-1= ρi-1/ρi-2p(i)= z(i-1) + βi-1 p(i-1)
endifq(i)= [A]p(i)αi = ρi-1/p(i)q(i)x(i)= x(i-1) + αip(i)r(i)= r(i-1) - αiq(i)check convergence |r|end

Ax − bDNorm2 r
Resid == =≤ Tol
BNorm2 bb