有限要素法プログラミング演習積分の計算

有限要素法プログラミング演習 積分の計算
2009 年 12 月 14 日
1
はじめに
• この演習では 2 次元トラス、三角形要素、四角形アイソパラメトリック要素
を取り扱うが、どの要素でも Cauchy の第一運動法則から導かれる σij δεij と
いう被積分関数を要素内で積分するという手続きが必要である。トラス要素
や、３節点三角形要素ではこの被積分関数が一定値になるので積分の外に出
すことができ、領域積分するといっても、要素の体積すなわちトラスであれ
ば、断面積 × 長さ、三角形であれば、面積 × 板厚をかけるだけなので、た
いした手間ではないが、アイソパラメトリック要素をプログラミングすると
きには、かなりの比重を占めることになる。
実際に要素剛性マトリックスの作成に入る前に、要素の種別による領域積分
の違いを理解するために円の面積や球の体積といった簡単な例題を有限要素
法と同じ手順でプログラムする。
2
トラス要素の場合
• 今回のトラス要素は要素内では断面積は一定、すなわち、断面形状が一定で
あることを仮定しているので、球であれば断面が円形のトラスを用いること
にする。円の面積を有限要素法で求めるといっても、断面の半径を r とし
て、解析的に πr 2 で求めることになる。
• これを元に球の体積を求めてみよう。構造物の解析を行う、という観点から
は現実的ではないが、例えば、球を図に示すような、薄い円盤のスライスを
m 枚積み重ねたものとしてモデル化する。
区間 [0, r] を m 等分して、xi = r/m × (i − 1), x1 = 0, xm+1 = r とする。こ
のとき、i 枚目の円盤の半径を xi とすると、円盤の厚さは r/m 円盤の半径
は r 2 − x2i なので、球の体積 V は、以下のように近似できる。
m
2 r
r
2
·π r −
(i − 1)
V ≈
m
m
i=1
(1)
• このコーディング例を truss_int.f, truss_int.c とする。プログラム中
では体積は vol, 半径 r, 分割数 mp, 円周率 pi としている。
!########################################################################
!
truss_int.f
!########################################################################
implicit none
integer :: mp, i
real(8) :: pi, vol, r
r
mp
pi
vol
=
=
=
=
2.D0
300
3.14159265358979323846264338327950288
0.D0
do i = 1, mp
3
vol = vol + r / mp * pi* (r**2 - (r / mp * (i-1))**2 )
end do
write(*,*) vol*2.D0, 4.D0 * pi * r**3 / 3
stop
end
/************************************************************************
truss_int.c
************************************************************************/
#include <stdio.h>
main()
{
double pi, vol, r;
int mp, i;
r = 2.0;
mp = 300;
pi = 3.14159265358979323846264338327950288;
vol = 0.0;
for (i=0; i<mp; i++){
vol += r / mp * pi* (r*r - (r / mp * i)*(r / mp * i) );
}
printf("%lg, %lg\n", vol*2.0, 4.0*pi*r*r*r/3.0);
}
4
円の領域分割–メッシュ生成
• 次に三角形要素の場合を考える。すなわち、「三角形要素を使って球の体積
を求めよ。」という問題であるが、最初から球の体積だとイメージに飛躍が
ありすぎるので、まずは円の面積を求めてみることにする。といっても、や
はりイメージがわかない人が多いであろう。ここでは有限要素法の手順にし
たがって、順を追って面積を求めてみる。
円の面積 A を求める式を考えると、中心の座標 (0, 0) 半径 r の円形の領域
を V として、以下のように表すことができる。
A=
dX1 dX2
(2)
V
きわめて自明の式であるが、この積分領域を重なり合いのない小領域に分割
して、個々の領域で積分を実行して、後から和を取る、という操作を行うと、
以下のように表すことができる。
A=
dX1 dX2 =
dX1 dX2
(3)
V
e
Ve
ここで領域を図のように三角形の集まりとして分割した場合を考える。結果
的には、円の外周を折れ線で近似することになるので、厳密な面積を求める
ことができないが、それでも分割数を細かくすると、分割数を細かくしてい
けばよい近似値が得られると期待できそうである。
• ところが、この三角形に分割するプログラムを作成せよだけでもやはり大変
になりそうな予感がするだろう。実際にこの三角形に分割するという作業は
一般にはメッシュ生成と呼ばれ、設計実務においても完全な自動化は困難で、
メッシュ生成だけでひとつの研究分野になっているぐらいである。ただ、簡
単な例題に限ればデータの生成だけでもよい例題になる。ここでは、後述す
5
る四角形要素でも同じメッシュ生成プログラムを使うことにして、まずは領
域を四角形に分割して、それぞれの四角形を２分割することにより三角形分
割を求めよう。なお、このプログラム作成は、本筋から若干離れているので、
完成版のプログラムを利用できるだけでもよい。
6
節点の生成
• このような領域のメッシュ分割には定石がある。四角形領域であればたとえ正
方形でなくとも各辺を均等に分割して、対辺とつなげばある程度形のそろっ
た四角形に分割できる。形がそろっていることは有限要素法で解析を行った
ときの精度の良し悪しに反映されるので、メッシュを生成する際には、十分
注意する。この四角形を半分に分割すれば三角形要素によるメッシュも生成
できる。
三角形の領域の場合は、重心と各辺の中点を結ぶと四角形に分割することが
でき、この四角形を前述のように細分して、出来上がった四角形を半分に分
割することにより、三角形要素によるメッシュが生成できる。
• 円形領域といっても、基本的には三角形の変形したものととらえる。まず、
図のように P1 , · · · P7 をとる。
四角形 P1 − P2 − P4 − P6 内部に均等に節点を配置する。x1 , x2 方向同時に
ループインデックスで表現するとバグの温床になるので、P1 − P2 , P6 − P4 間
に中間節点 Q1 , Q2 を設置する。ただし i = 0, · · · , m であることに注意する。
i
(P2 − P1 )
m
i
Q2 = P6 + (P4 − P6 )
m
Q1 = P1 +
7
(4)
(5)
P7 , (0, r)
P5 ,
P6 , 0, r2 , 0
P4 ,
P1 , (0, 0)
r
√
, √r
2 2 2 2
P2 ,
r
2
,0
√r , √r
2
2
P3 , (r, 0)
求める節点座標 Cij は同様に j = 0, · · · , m として以下のように計算できる。
Cij = Q1 +
j
(Q2 − Q1 )
m
P6 , 0, 2r , 0
Q2
P4 ,
r
√
, √r
2 2 2 2
Cij
P1 , (0, 0)Q1
P3
8
(6)
• 次に扇形 P2 −P3 −P5 −P7 −P6 −P4 の内部に放射状に節点を配置する。P2 −P3
を m 等分して同心円の半径を下式によって求める。ただし i = 1, · · · , m で
ある。
i
S = P2 + (P3 − P2 )
(7)
m
角度方向は P2 から P6 に向かう方向を正にとり、j = 0, · · · 2m として以下
のように与えられる。
π
j
(8)
θ=
2(2m)
これらを用いると求める節点座標 Cij は
C1 = S cos θ
(9)
C2 = S sin θ
(10)
P7
P5
P6
P4
S
C
θ
P1
P2
9
P3
可視化データの作成 基礎編
• この結果が正しいかどうかを確認するのには、可視化すると直感的に理解で
きる。ここでは MicroAVS を使った可視化のためのデータの作成について説
明する。MicroAVS はさまざまな形式のデータを可視化できるが、有限要素
法のデータは「非構造格子型データ」に分類される。拡張子は .inp。 形の
そろっていないポリゴンの集合と考えればイメージできる人も多いと思う。
メッシュデータを表示するだけなら１ステップ分のデータで十分であるが、
ソフトは複数ステップのデータに対応していて物体の動作なども可視化でき
るようになっている。
フォーマットは以下のようになっている。これだけでも機能が多いので、必
要な箇所だけ説明する。詳しくはソフトに付属のマニュアルが親切なので参
照されたい。
( コメント行 )
( ステップ数 )
( データの繰り返しタイプ )
( ステップ番号 1 ) ( コメント )
( 全節点数 ) ( 全要素数 )
( 節点番号 1 ) ( X 座標 ) ( Y 座標 ) ( Z 座標 )
( 節点番号 2 ) ( X 座標 ) ( Y 座標 ) ( Z 座標 )
・
・
( 要素番号 1 ) ( 材料番号 ) ( 要素の種類 ) ( 要素を構成する節点のつながり )
( 要素番号 2 ) ( 材料番号 ) ( 要素の種類 ) ( 要素を構成する節点のつながり )
・
・
( 節点のデータ数 ) ( 要素のデータ数 )
・
・
( 節点のデータ成分数 ) ( 節点成分 1 の構成数 ) ( 節点成分 2 の構成数 ) ・
( 節点データ成分 1 のラベル )，( 単位 )
( 節点データ成分 2 のラベル )，( 単位 )
・
・
( 各節点データ成分のラベル )，( 単位 )
・
・
・
・
( 節点番号 1 ) ( 節点データ 1 ) ( 節点データ 2 ) ・
( 節点番号 2 ) ( 節点データ 1 ) ( 節点データ 2 ) ・
・
・
・
・
・
・
( 要素のデータ成分数 ) ( 要素成分 1 の構成数 ) ( 要素成分 2 の構成数 ) ・
・
・
10
( 要素データ成分 1 のラベル )，( 単位 )
( 要素データ成分 2 のラベル )，( 単位 )
・
・
( 要素番号 1 ) ( 要素データ 1 ) ( 要素データ 2 ) ・
・
・
・
・
( 要素番号 2 ) ( 要素データ 1 ) ( 要素データ 2 ) ・
・
・
・
・
・
・
( ステップ番号 2 ) ( コメント )
( 全節点数 ) (全要素数 )
以下繰り返し
• コメント行:今回は必要なし、かつ省略可
• ステップ数:本来はアニメーションのコマ送りのように各ステップで 1 画像を
作成し連続した動きを表示するためのものだが、今回はメッシュを表示する
だけなので 1 でよい。
• データの繰り返しタイプ:複数データを表示するときに読み替えるデータ。形
状が変化せずデータだけ変化していく場合は data. 形状のみ変化する場合
は geom 両方とも変化する場合は data_geom のキーワードを指定する。１ス
テップしかない場合は data_geom を指定する。
• 全節点数、全要素数:これは可視化する節点数と要素数を指定する。今回は節
点しかないが、要素タイプが「点要素」であるとして要素数も節点数と同じ
数を指定する。
• 節点番号、X 座標、Y 座標、Z 座標:3 次元データのみに対応しているので、
2 次元データは Z 座標を 0 として指定する。
• 要素番号、材料番号、要素の種類、要素を構成する節点のつながり:今回は
材料番号は使用しないので今後も含めて 1 を指定すればよい。要素の種類
は、トラスなら line 三角形なら tri 四角形なら quad 点のデータならば pt
を指定する。要素を構成する節点のつながりは、講義で用いているものと同
じで、三角形、四角形とも反時計回りに指定する。トラスはどちらを先に記
述しても可視化のみなので影響はない。点の場合は節点番号を指定する。
• 節点のデータ数、要素のデータ数:今回の演習では 1 0 を指定する。
• 節点のデータ成分数、節点成分 1 の構成数、節点成分 2 の構成数 ・
・
・:上記
のように節点のデータ数、要素のデータ数を 1 0 と指定した場合は、1 1 を
指定する。
11
• 節点データ成分 1 のラベル，単位:今回は”,” のみを指定する。
• 節点番号、節点データ 1、節点データ 2 ・
・
・:今回は節点データとして節点番
号を記入する。
• 要素のデータ成分数以降のデータ:今回は必要なし。省略可
12
節点生成のサンプルプログラム
• 例えば分割数 m = 3 として、この段階までのデータを生成するプログラム
を mesh1.f, mesh1.c とする。また、可視化データのファイルを point.inp（ダ
ウンロードのみ） とする。
37
36
35
30
29
28
23
34
22
27
21
4
3
2
1
8
12
6
5
10
9
26
19
16
11 15
7
33
20
32
18
14
13
17
25
24
31
!######################################################################
!
mesh1.f
!######################################################################
implicit none
real(8) :: coords(2,10000), p1(2), p2(2), p3(2), p4(2), p5(2),
$
p6(2), p7(2), q1(2), q2(2), r, pi,s
integer :: i, j, mp, nnode, nelem, lnods(4,10000)
pi = 3.14159265358979323846264338327950288
r = 2.D0
mp = 3
p1(1)
p1(2)
p2(1)
p2(2)
p3(1)
p3(2)
p6(1)
p6(2)
p7(1)
p7(2)
=
=
=
=
=
=
=
=
=
=
0.D0
0.D0
r / 2.D0
0.D0
r
0.D0
0.D0
r / 2.D0
0.D0
r
p4(1) = dsqrt(2.D0) / 4.D0 * r
p4(2) = dsqrt(2.D0) / 4.D0 * r
nnode = 0
do i = 0, mp
13
q1(1) = p1(1) + (p2(1)-p1(1)) / mp * i
q1(2) = p1(2) + (p2(2)-p1(2)) / mp * i
q2(1) = p6(1) + (p4(1)-p6(1)) / mp * i
q2(2) = p6(2) + (p4(2)-p6(2)) / mp * i
do j = 0, mp
nnode = nnode + 1
coords(1,nnode) = q1(1) + (q2(1)-q1(1)) / mp * j
coords(2,nnode) = q1(2) + (q2(2)-q1(2)) / mp * j
end do
end do
do i = 1, mp
s = p2(1) + (p3(1)
do j = 0, 2*mp
nnode = nnode +
coords(1,nnode)
coords(2,nnode)
end do
end do
- p2(1)) / mp * i
1
= s * dcos(pi / 2.D0 / (2*mp) * j)
= s * dsin(pi / 2.D0 / (2*mp) * j)
! debug write -----------------------------------------------------------open(10,file=’point.inp’)
write(10,*) 1
write(10,*) ’data_geom’
write(10,*) ’step1’
write(10,*) nnode, nnode
do i = 1, nnode
write(10,1001) i, coords(1:2,i), 0.D0
end do
1001 format(i4,x,e13.6,x,e13.6,x,e13.6)
do i = 1, nnode
write(10,1002) i, i
end do
1002 format(i4,’ 1 pt’,4i8)
write(10,*) 1, 0
write(10,*) 1, 1
write(10,*) ’,’
do i = 1, nnode
write(10,1001) i, dble(i)
end do
close(10)
stop
end
/**********************************************************************
mesh1.c
**********************************************************************/
#include <stdio.h>
main()
{
double coords[2][10000], p1[2], p2[2], p3[2], p4[2], p5[2],
p6[2], p7[2], q1[2], q2[2], r, pi, s;
int i, j, mp, nnode, nelem, lnods[4][10000];
FILE *fp;
pi = 3.14159265358979323846264338327950288;
r = 2.0;
mp = 3;
14
p1[0]
p1[1]
p2[0]
p2[1]
p3[0]
p3[1]
p6[0]
p6[1]
p7[0]
p7[1]
=
=
=
=
=
=
=
=
=
=
0.0;
0.0;
r / 2.0;
0.0;
r;
0.0;
0.0;
r / 2.0;
0.0;
r;
p4[0] = sqrt(2.0) / 4.0 * r;
p4[1] = sqrt(2.0) / 4.0 * r;
nnode = 0;
for(i=0; i<=mp; i++){
q1[0] = p1[0] + (p2[0]-p1[0]) / mp * i;
q1[1] = p1[1] + (p2[1]-p1[1]) / mp * i;
q2[0] = p6[0] + (p4[0]-p6[0]) / mp * i;
q2[1] = p6[1] + (p4[1]-p6[1]) / mp * i;
for(j=0; j<= mp; j++){
coords[0][nnode] = q1[0] + (q2[0]-q1[0]) / mp * j;
coords[1][nnode] = q1[1] + (q2[1]-q1[1]) / mp * j;
nnode = nnode + 1;
}
}
for(i=0; i<mp; i++){
s = p2[0] + (p3[0] - p2[0]) / mp * (i+1);
for(j=0; j<=mp*2; j++){
coords[0][nnode] = s * cos(pi / 2.0 / (2*mp) * j);
coords[1][nnode] = s * sin(pi / 2.0 / (2*mp) * j);
nnode = nnode + 1;
}
}
/* debug write -------------------------------------------------- */
if ((fp = fopen("point.inp","w")) == NULL)
exit(1);
fprintf(fp,
fprintf(fp,
fprintf(fp,
fprintf(fp,
"1\n");
"data_geom\n");
"step1\n");
"%d %d\n", nnode, nnode);
for(i=0; i<nnode ; i++){
fprintf(fp, "%4d %14.5e %14.5e %14.5e\n", i+1,
coords[0][i], coords[1][i], 0.0);
}
for(i=0; i<nnode ; i++){
fprintf(fp, "%4d 1 pt %8d\n", i,i);
}
fprintf(fp, "1 0\n");
fprintf(fp, "1 1\n");
fprintf(fp, ",\n");
for(i=0; i<nnode ; i++){
fprintf(fp, "%4d %14.5e\n", i+1, (double)(i+1) );
}
}
15
コネクティビティの生成
• 要素を構成する節点のつながり（コネクティビティ）を作成する。配列名は
lnods を使う。まず最初は考えやすい部分から作成して、最後に個別の処理
を必要とする箇所を考える。原点を含む要素のコネクティビティは (1, 4, 5, 2)
であるが、これは分割数 m を用いて表せば、(1, 1 + (m + 1), 2 + (m + 1), 2)
である。これを要素 1 として、上の要素を 2 とすると、要素 2 のコネクティ
ビティは要素 1 のコネクティビティにすべて 1 加えたものであることがわか
る。これより、上方向のループカウンターを j とすると、j = 0, · · · , m − 1
について
lnods(1,nelem) = 1
+ j
lnods(2,nelem) = 1 + (m+1) + j
lnods(3,nelem) = 2 + (m+1) + j
lnods(4,nelem) = 2
+ j
となる。要素 1 の右側の要素のコネクティビティは (5, 9, 10, 6) であるが、こ
れは要素 1 のコネクティビティにすべて (m + 1) 加えたものであることがわ
かる。(1 + (m + 1), 4 + (m + 1), 5 + (m + 1), 2 + (m + 1)) ループカウンター
を i とすると、i = 0, · · · , m − 1 について
lnods(1,nelem) = 1
+ j + (m+1)*i
lnods(2,nelem) = 1 + (m+1) + j + (m+1)*i
lnods(3,nelem) = 2 + (m+1) + j + (m+1)*i
lnods(4,nelem) = 2
+ j + (m+1)*i
となる。
次に、外周部分の (17, 24, 25, 18) で構成される要素から考える。これを分割
数 m を用いて表すと、以下のようになる。
lnods(1,nelem) = (m+1)*(m+1) + 1
lnods(2,nelem) = (m+1)*(m+1) + 1 + (2*m+1)
lnods(3,nelem) = (m+1)*(m+1) + 2 + (2*m+1)
lnods(4,nelem) = (m+1)*(m+1) + 2
したがって、周方向のループカウンターを j 半径方向のループカウンターを
i とすると、i = 0, · · · , m − 2, j = 0, · · · , 2m − 1 について
lnods(1,nelem) = (m+1)*(m+1) + 1
+ j + (2*m+1)*i
lnods(2,nelem) = (m+1)*(m+1) + 1 + (2*m+1) + j + (2*m+1)*i
lnods(3,nelem) = (m+1)*(m+1) + 2 + (2*m+1) + j + (2*m+1)*i
lnods(4,nelem) = (m+1)*(m+1) + 2
+ j + (2*m+1)*i
最後に、中央の四角形部分と、外周部分をつなぐ 1 層分のコネクティビティ
を考える。(13, 17, 18, 14) は分割数 m を用いて表すと、以下のようになる。
16
lnods(1,nelem) = (m+1)*m + 1
lnods(2,nelem) = (m+1)*(m+1) + 1
lnods(3,nelem) = (m+1)*(m+1) + 2
lnods(4,nelem) = (m+1)*m + 1
この図の例であれば 3(= m) 要素分が同じ順に節点が並んでいる。ループカウ
ンターを i とし、て i = 0, m−1 について lnods(1,nelem) = (m+1)*m + 1
+ i
lnods(2,nelem) = (m+1)*(m+1) + 1 + i
lnods(3,nelem) = (m+1)*(m+1) + 2 + i
lnods(4,nelem) = (m+1)*m + 1
+ i
(16, 20, 21, 12) の要素からは、パターンが変わるが、まずこの要素のコネク
ティビティを分割数 m を使って表すと、
lnods(1,nelem) = (m+1)*(m+1) lnods(2,nelem) = (m+1)*(m+1) + m+1
lnods(3,nelem) = (m+1)*(m+1) + m+2 lnods(4,nelem) = (m+1)* m
左
隣に 3(= m) 要素分が同じ順に節点が並んでいる。ループカウンターを i とし
て i = 0, m−1 について lnods(1,nelem) = (m+1)*(m+1)
- (m+1)*i
lnods(2,nelem) = (m+1)*(m+1) + m+1 + i lnods(3,nelem) = (m+1)*(m+1) + m+2 + i
lnods(4,nelem) = (m+1)* m
- (m+1)*i
37
36
35
30
29
28
23
34
22
27
21
4
33
8
20
12
3
7
11
2
1
6
5
26
19
16
15
32
18
10
9
14
13
17
25
24
31
• 四角形要素のコネクティビティを分割して、三角形要素のコネクティビティを
17
作成する。1 つの四角形要素を 2 つの三角形要素に分割する。出来上がった
コネクティビティの配列 lnods をもとに、三角形のコネクティビティ lnods3
を構成する。具体的には、
lnods3(1,nelem3) = lnods(1,nelem)
lnods3(2,nelem3) = lnods(2,nelem)
lnods3(3,nelem3) = lnods(4,nelem)
nelem3 ++
lnods3(1,nelem3) = lnods(2,nelem)
lnods3(2,nelem3) = lnods(3,nelem)
lnods3(3,nelem3) = lnods(4,nelem)
• またこれを可視化するための MicroAVS のデータを作成する。それぞれ quad.inp,
tri.inp（ダウンロードのみ） とする。先に説明したデータ形式で、要素のデー
タを記述する部分が置き換わるだけで、後は変更しなくても大丈夫である。
具体的には、要素の種類を、pt から quad, tri に変更し、それぞれ４節点、
３節点分の節点番号を記述する。
• 最後に、データが正しくできていることを確認できたら、三角形と、四角形
のデータを tri.dat, quad.dat （ダウンロードのみ）として書き出しておく。
書き出しフォーマットは、対応する読み込み部分を後からプログラムするの
で、独自に設定してかまわない。もちろん MicroAVS のデータ形式でもかま
わないが、サンプルプログラムでは、必要最低限だけを書き出すようにして
ある。
18
37
36
35
30
29
28
23
34
22
27
21
4
33
8
20
12
3
7
11
2
6
1
26
19
16
15
32
18
10
5
9
14
13
17
25
24
31
メッシュ生成のサンプルプログラム
• これまで説明してきたことをまとめたサンプルコードを示す。
!########################################################################
!
mesh.f
!########################################################################
implicit none
real(8) :: coords(2,10000), p1(2), p2(2), p3(2), p4(2), p5(2),
$
p6(2), p7(2), q1(2), q2(2), r, pi, s
integer :: i, j, mp, nnode, nelem, lnods(4,10000),
$
lnods3(3,10000), nelem3
pi = 3.14159265358979323846264338327950288
r = 2.D0
mp = 3
p1(1)
p1(2)
p2(1)
p2(2)
p3(1)
p3(2)
p6(1)
p6(2)
p7(1)
p7(2)
=
=
=
=
=
=
=
=
=
=
0.D0
0.D0
r / 2.D0
0.D0
r
0.D0
0.D0
r / 2.D0
0.D0
r
19
p4(1) = dsqrt(2.D0) / 4.D0 * r
p4(2) = dsqrt(2.D0) / 4.D0 * r
nnode = 0
do i = 0, mp
q1(1) = p1(1) + (p2(1)-p1(1)) / mp * i
q1(2) = p1(2) + (p2(2)-p1(2)) / mp * i
q2(1) = p6(1) + (p4(1)-p6(1)) / mp * i
q2(2) = p6(2) + (p4(2)-p6(2)) / mp * i
do j = 0, mp
nnode = nnode + 1
coords(1,nnode) = q1(1) + (q2(1)-q1(1)) / mp * j
coords(2,nnode) = q1(2) + (q2(2)-q1(2)) / mp * j
end do
end do
do i = 1, mp
s = p2(1) + (p3(1)
do j = 0, 2*mp
nnode = nnode +
coords(1,nnode)
coords(2,nnode)
end do
end do
- p2(1)) / mp * i
1
= s * dcos(pi / 2.D0 / (2*mp) * j)
= s * dsin(pi / 2.D0 / (2*mp) * j)
! debug write -----------------------------------------------------------open(10,file=’point.inp’)
write(10,*) 1
write(10,*) ’data_geom’
write(10,*) ’step1’
write(10,*) nnode, nnode
do i = 1, nnode
write(10,1001) i, coords(1:2,i), 0.D0
end do
1001 format(i4,x,e13.6,x,e13.6,x,e13.6)
do i = 1, nelem
write(10,1002) i, i
end do
1002 format(i4,’ 1 pt’,4i8)
write(10,*) 1, 0
write(10,*) 1, 1
write(10,*) ’,’
do i = 1, nnode
write(10,1000) i, 0.D0
end do
close(10)
! connectivity -----------------------------------------------------------nelem = 0
do i = 0, mp-1
do j = 0, mp-1
nelem = nelem
lnods(1,nelem)
lnods(2,nelem)
lnods(3,nelem)
lnods(4,nelem)
end do
end do
+
=
=
=
=
1
1
+ j + i * (mp+1)
1 + (mp+1) + j + i * (mp+1)
2 + (mp+1) + j + i * (mp+1)
2
+ j + i * (mp+1)
20
do i = 0, mp-2
do j = 0, 2*mp-1
nelem = nelem + 1
lnods(1,nelem) = (mp+1)**2
lnods(2,nelem) = (mp+1)**2
lnods(3,nelem) = (mp+1)**2
lnods(4,nelem) = (mp+1)**2
end do
end do
do i = 0, mp-1
nelem = nelem + 1
lnods(1,nelem) = (mp+1)*mp
lnods(2,nelem) = (mp+1)**2
lnods(3,nelem) = (mp+1)**2
lnods(4,nelem) = (mp+1)*mp
end do
+
+
+
+
+
+
+
+
1
1
2
2
+
+
+
+
1
+
1 + (2*mp+1) +
2 + (2*mp+1) +
2
+
j
j
j
j
+
+
+
+
i
i
i
i
i
i
i
i
do i = 0, mp-1
nelem = nelem + 1
lnods(1,nelem) = (mp+1)**2
lnods(2,nelem) = (mp+1)**2 + mp+1 +
lnods(3,nelem) = (mp+1)**2 + mp+2 +
lnods(4,nelem) = (mp+1)*mp
end do
(mp+1)*i
i
i
(mp+1)*i
nelem3 = 0
do i = 1, nelem
nelem3 = nelem3 + 1
lnods3(1,nelem3) = lnods(1,i)
lnods3(2,nelem3) = lnods(2,i)
lnods3(3,nelem3) = lnods(4,i)
nelem3 = nelem3 + 1
lnods3(1,nelem3) = lnods(2,i)
lnods3(2,nelem3) = lnods(3,i)
lnods3(3,nelem3) = lnods(4,i)
end do
! debug write quad element -----------------------------open(10,file=’quad.inp’)
write(10,*) 1
write(10,*) ’data_geom’
write(10,*) ’step1’
write(10,*) nnode, nelem
do i = 1, nnode
write(10,1000) i, coords(1:2,i), 0.D0
end do
1000 format(i4,x,e13.6,x,e13.6,x,e13.6)
do i = 1, nelem
write(10,2000) i, lnods(1:4,i)
end do
2000 format(i4,’ 1 quad’,4i8)
write(10,*) 1, 0
write(10,*) 1, 1
write(10,*) ’,’
do i = 1, nnode
write(10,1000) i, 0.D0
end do
close(10)
! debug write triangel element -----------------------------open(10,file=’tri.inp’)
21
*
*
*
*
(2*mp+1)
(2*mp+1)
(2*mp+1)
(2*mp+1)
write(10,*) 1
write(10,*) ’data_geom’
write(10,*) ’step1’
write(10,*) nnode, nelem*2
do i = 1, nnode
write(10,1000) i, coords(1:2,i), 0.D0
end do
do i = 1, nelem3
write(10,2001) i, lnods3(1:3,i)
end do
2001 format(i4,’ 1 tri’,4i8)
write(10,*) 1, 0
write(10,*) 1, 1
write(10,*) ’,’
do i = 1, nnode
write(10,1000) i, 0.D0
end do
close(10)
! data file triangle -------------------------------------------------open(20,file=’tri.dat’)
write(20,*) nnode, nelem3
do i = 1, nnode
write(20,1000) i, coords(1:2,i)
end do
do i = 1, nelem3
write(20,3000) i, lnods3(1:3,i)
end do
close(20)
! data file quad -------------------------------------------------open(20,file=’quad.dat’)
write(20,*) nnode, nelem
do i = 1, nnode
write(20,1000) i, coords(1:2,i)
end do
do i = 1, nelem
write(20,3000) i, lnods(1:4,i)
end do
3000 format(5i8)
close(20)
stop
end
/**********************************************************************
mesh.c
**********************************************************************/
#include <stdio.h>
main()
{
double coords[2][10000], p1[2], p2[2], p3[2], p4[2], p5[2],
p6[2], p7[2], q1[2], q2[2], r, pi, s;
int i, j, mp, nnode, nelem, lnods[4][10000], lnods3[3][10000], nelem3;
FILE *fp;
pi = 3.14159265358979323846264338327950288;
r = 2.0;
22
mp = 3;
p1[0]
p1[1]
p2[0]
p2[1]
p3[0]
p3[1]
p6[0]
p6[1]
p7[0]
p7[1]
=
=
=
=
=
=
=
=
=
=
0.0;
0.0;
r / 2.0;
0.0;
r;
0.0;
0.0;
r / 2.0;
0.0;
r;
p4[0] = sqrt(2.0) / 4.0 * r;
p4[1] = sqrt(2.0) / 4.0 * r;
nnode = 0;
for(i=0; i<=mp; i++){
q1[0] = p1[0] + (p2[0]-p1[0]) / mp * i;
q1[1] = p1[1] + (p2[1]-p1[1]) / mp * i;
q2[0] = p6[0] + (p4[0]-p6[0]) / mp * i;
q2[1] = p6[1] + (p4[1]-p6[1]) / mp * i;
for(j=0; j<= mp; j++){
coords[0][nnode] = q1[0] + (q2[0]-q1[0]) / mp * j;
coords[1][nnode] = q1[1] + (q2[1]-q1[1]) / mp * j;
nnode = nnode + 1;
}
}
for(i=0; i<mp; i++){
s = p2[0] + (p3[0] - p2[0]) / mp * (i+1);
for(j=0; j<=mp*2; j++){
coords[0][nnode] = s * cos(pi / 2.0 / (2*mp) * j);
coords[1][nnode] = s * sin(pi / 2.0 / (2*mp) * j);
nnode = nnode + 1;
}
}
/* debug write -------------------------------------------------- */
if ((fp = fopen("point.inp","w")) == NULL)
exit(1);
fprintf(fp,
fprintf(fp,
fprintf(fp,
fprintf(fp,
"1\n");
"data_geom\n");
"step1\n");
"%d %d\n", nnode, nnode);
for(i=0; i<nnode ; i++){
fprintf(fp, "%4d %14.5e %14.5e %14.5e\n", i+1,
coords[0][i], coords[1][i], 0.0);
}
for(i=0; i<nnode ; i++){
fprintf(fp, "%4d 1 pt %8d\n", i,i);
}
fprintf(fp, "1 0\n");
fprintf(fp, "1 1\n");
fprintf(fp, ",\n");
for(i=0; i<nnode ; i++){
fprintf(fp, "%4d %14.5e\n", i+1, (double)(i+1) );
}
/* connectivity ------------------------------------------------------------*/
23
nelem = 0;
for (i=0; i<=mp-1; i++){
for(j=0; j<=mp-1; j++){
lnods[0][nelem] = 1
+ j
lnods[1][nelem] = 1 + (mp+1) + j
lnods[2][nelem] = 2 + (mp+1) + j
lnods[3][nelem] = 2
+ j
nelem ++;
}
}
+
+
+
+
i
i
i
i
*
*
*
*
(mp+1);
(mp+1);
(mp+1);
(mp+1);
for (i=0; i<=mp-2; i++){
for (j=0; j<=2*mp-1; j++){
lnods[0][nelem] = (mp+1)*(mp+1)
lnods[1][nelem] = (mp+1)*(mp+1)
lnods[2][nelem] = (mp+1)*(mp+1)
lnods[3][nelem] = (mp+1)*(mp+1)
nelem ++;
}
}
+
+
+
+
1
+
1 + (2*mp+1) +
2 + (2*mp+1) +
2
+
for (i=0; i<=mp-1; i++){
lnods[0][nelem] = (mp+1)*mp
lnods[1][nelem] = (mp+1)*(mp+1)
lnods[2][nelem] = (mp+1)*(mp+1)
lnods[3][nelem] = (mp+1)*mp
nelem ++;
}
1
1
2
2
+
+
+
+
+
+
+
+
j
j
j
j
+
+
+
+
i
i
i
i
*
*
*
*
i;
i;
i;
i;
for (i=0; i<=mp-1; i++){
lnods[0][nelem] = (mp+1)*(mp+1)
lnods[1][nelem] = (mp+1)*(mp+1) + mp+1 +
lnods[2][nelem] = (mp+1)*(mp+1) + mp+2 +
lnods[3][nelem] = (mp+1)*mp
nelem ++;
}
(mp+1)*i;
i;
i;
(mp+1)*i;
nelem3 = 0;
for(i=0; i<nelem; i++){
lnods3[0][nelem3] = lnods[0][i];
lnods3[1][nelem3] = lnods[1][i];
lnods3[2][nelem3] = lnods[3][i];
nelem3 ++;
lnods3[0][nelem3] = lnods[1][i];
lnods3[1][nelem3] = lnods[2][i];
lnods3[2][nelem3] = lnods[3][i];
nelem3 ++;
}
/* debug write quad element ------------------------------*/
if ((fp = fopen("quad.inp","w")) == NULL)
exit(1);
fprintf(fp,
fprintf(fp,
fprintf(fp,
fprintf(fp,
"1\n");
"data_geom\n");
"step1\n");
"%d %d\n", nnode, nelem);
for(i=0; i<nnode ; i++){
fprintf(fp, "%4d %14.5e %14.5e %14.5e\n", i+1,
coords[0][i], coords[1][i], 0.0);
}
for(i=0; i<nelem ; i++){
fprintf(fp, "%4d 1 quad %8d %8d %8d %8d\n", i+1,
lnods[0][i],lnods[1][i],lnods[2][i],lnods[3][i]);
24
(2*mp+1);
(2*mp+1);
(2*mp+1);
(2*mp+1);
}
fprintf(fp, "1 0\n");
fprintf(fp, "1 1\n");
fprintf(fp, ",\n");
for(i=0; i<nnode ; i++){
fprintf(fp, "%4d %14.5e\n", i+1, (double)(i+1) );
}
fclose(fp);
/* debug write triangel element ------------------------------*/
if ((fp = fopen("tri.inp","w")) == NULL)
exit(1);
fprintf(fp,
fprintf(fp,
fprintf(fp,
fprintf(fp,
"1\n");
"data_geom\n");
"step1\n");
"%d %d\n", nnode, nelem3);
for(i=0; i<nnode ; i++){
fprintf(fp, "%4d %14.5e %14.5e %14.5e\n", i+1,
coords[0][i], coords[1][i], 0.0);
}
for(i=0; i<nelem3 ; i++){
fprintf(fp, "%4d 1 tri %8d %8d %8d\n", i+1,
lnods3[0][i],lnods3[1][i],lnods3[2][i]);
}
fprintf(fp, "1 0\n");
fprintf(fp, "1 1\n");
fprintf(fp, ",\n");
for(i=0; i<nnode ; i++){
fprintf(fp, "%4d %14.5e\n", i+1, (double)(i+1) );
}
fclose(fp);
/* data file triangle --------------------------------------------------*/
if ((fp = fopen("tri.dat","w")) == NULL)
exit(1);
fprintf(fp, "%d %d\n", nnode, nelem3);
for(i=0; i<nnode ; i++){
fprintf(fp, "%4d %14.5e %14.5e\n", i+1,
coords[0][i], coords[1][i]);
}
for(i=0; i<nelem3 ; i++){
fprintf(fp, "%4d %8d %8d %8d\n", i+1,
lnods3[0][i],lnods3[1][i],lnods3[2][i]);
}
fclose(fp);
/* data file quad --------------------------------------------------*/
if ((fp = fopen("quad.dat","w")) == NULL)
exit(1);
fprintf(fp, "%d %d\n", nnode, nelem);
for(i=0; i<nnode ; i++){
fprintf(fp, "%4d %14.5e %14.5e\n", i+1,
coords[0][i], coords[1][i]);
}
for(i=0; i<nelem ; i++){
fprintf(fp, "%4d %8d %8d %8d %8d\n", i+1,
25
lnods[0][i],lnods[1][i],lnods[2][i],lnods[3][i]);
}
fclose(fp);
}
26
三角形要素による円の面積
• 円のメッシュ生成ができたので、今度は面積を求めてみる。あたらためて、
有限要素法の手順に従って円の面積 A を求める式を考えると、中心の座標
(0, 0) 半径 r の円形の領域を V として、以下のように表すことができる。
A=
dX1 dX2
(11)
V
この積分領域を重なり合いのない三角形あるいは四角形の小領域に分割して、
個々の領域で積分を実行して、後から和を取る、という操作を行うと、以下
のように表すことができる。
A=
dX1 dX2 =
dX1 dX2
(12)
V
Ve
e
• 先に作成したメッシュを用いて実際にこれを行ってみる。メッシュ生成プロ
グラムの最後に面積を計算する部分を付け加えてもよいが、ファイルが長く
なると、それだけでやはりバグの温床になるので、、まずは、データファイ
ルから三角形の節点座標とコネクティビティを読み込むことにする。実は後
述する球の体積を求めるプログラムではメッシュ分割に伴って誤差が減少し
ていく様子を見るために、メッシュ生成と計算プログラムを合体することに
なる。
• 最初から複数要素のデータをファイルから読み込んで円の面積を計算すると、
やはりこれもバグの温床になるので、まずは確実に三角形の面積を求める。
(1)
(1)
(2)
(2)
(3)
(3)
三角形の面積は、三角形の節点座標を X1 , X2 , X1 , X2 , X1 , X2
て、
2A =
(2)
(1)
X 1 − X1
(2)
(1)
X 2 − X2
(2)
(3)
(1)
X1 − X1
×
(3)
(1)
X2 − X 2
(1)
(3)
(1)
(13)
(3)
(1)
(2)
(1)
= (X1 − X1 )(X2 − X2 ) − (X1 − X1 )(X2 − X2 )
(2)
(3)
(1)
(2)
(3)
(1)
とし
(2)
(1)
(3)
(2)
(14)
(1)
(3)
= X 1 X2 + X 1 X2 + X 1 X 2 − X 1 X2 − X 1 X2 − X 1 X2
(15)
• このままコーディングすると、インデックスが多すぎて逆にわかにりくく
(1)
(1)
(2)
(2)
なってしまう危険性があるので、X1 = x1 , X2 = y1 , X1 = x2 , X2 =
(3)
(3)
y2 , X1 = x3 , X2 = y3 と書き直すと、比較的読みやすいコードになる。
27
(3)
(3)
(3)(X1 , X2 )
(2)
(2)
(2)(X1 , X2 )
(1)
(2)
(1)(X1 , X2 )
2A =
x2 − x1
y2 − x1
x3 − x1
×
y3 − y1
(16)
= (x2 − x1 )(y3 − y1 ) − (x3 − x1 )(y2 − y1 )
(17)
= x2 y3 + x1 y2 + x3 y1 − x2 y1 − x3 y2 − x1 y3
(18)
• 例えば、(0, 0), (1, 0), (0, 1) や、これを平行移動した (1, 2), (2, 2), (1, 3) の面積
√
√
が正しく 12 と計算できることを確認したら、辺が軸に平行でない (0, 0), ( 3, 1), (1, 3)
(面積 1) や、(0, 0), (5, 2), (3, 5) (面積 9.5) についても確認すれば、面積の公
式のコーディングにバグが含まれてる可能性は低い。
• もしも値が異なっていたら、すなわち、バグがあったなら、以下のように２
通りの公式をコーディングして値を確認すると見つけやすい。
a = ((x2-x1)*(y3-y1)-(x3-x1)*(y2-y1)) / 2.0
a = (x2*y3 + x1*y2 + x3*y1 - x2*y1 - x3*y2 - x1*y3) / 2.0
数値計算は、とりわけ変数が多いので、変数名を定式化で用いた変数にあわ
せ、式も全く同じ式をそのままコーディングすると、比較的読みやすいコー
ドが作成でき、デバッグが容易になる。ただし、このようにコーディングす
ると、その引き換えに、コードが長くなってそのために読みにくくなるなど
の弊害があるため、好まない人も多いが、まずは正しく動作するコードを作
成した後に、コードの整理を行うようにするほうがデバッグに時間をとられ
ることが少ないようである。
サンプルコードを tri_area0.f, tri_area0.c とする。
!########################################################################
28
3
2
1
0
0
1
√
3
1
2
5
√
3−1
2
0
√
3−1
√
0
3
29
3
5
!
tri_area0.f
!########################################################################
implicit none
real(8) :: a, x1, y1, x2, y2, x3, y3
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
!
ex 1 -------------------------------------------------x1 = 0.D0
y1 = 0.D0
x2 = 1.D0
y2 = 0.D0
x3 = 0.D0
y3 = 1.D0
ex 2 -------------------------------------------------x1 = 0.D0 + 1.D0
y1 = 0.D0 + 2.D0
x2 = 1.D0 + 1.D0
y2 = 0.D0 + 2.D0
x3 = 0.D0 + 1.D0
y3 = 1.D0 + 2.D0
ex 3 -------------------------------------------------x1 = 0.D0
y1 = 0.D0
x2 = dsqrt(3.D0)
y2 = 1.D0
x3 = 1.D0
y3 = dsqrt(3.D0)
ex 4 -------------------------------------------------x1 = 0.D0
y1 = 0.D0
x2 = 5.D0
y2 = 2.D0
x3 = 3.D0
y3 = 5.D0
write(*,*) ’x1, y1’, x1, y1
write(*,*) ’x2, y2’, x2, y2
write(*,*) ’x3, y3’, x3, y3
!
a = ((x2-x1)*(y3-y1)-(x3-x1)*(y2-y1)) / 2.D0
a = (x2*y3 + x1*y2 + x3*y1 - x2*y1 - x3*y2 - x1*y3) / 2.D0
write(*,*) ’area’, a
stop
end
/**********************************************************************
tri_area0.c
**********************************************************************/
#include <stdio.h>
main()
{
double a, x1, y1, x2, y2, x3, y3;
/* ex1 -------------------------------------------------- */
/* x1 = 0.0;
y1 = 0.0;
x2 = 1.0;
y2 = 0.0;
x3 = 0.0;
y3 = 1.0; */
/*
/*
y1
x2
y2
ex 2
x1 =
= 0.0
= 1.0
= 0.0
-------------------------------------------------- */
0.0 + 1.0;
+ 2.0;
+ 1.0;
+ 2.0;
30
x3 = 0.0 + 1.0;
y3 = 1.0 + 2.0; */
/*
/*
y1
x2
y2
x3
y3
/*
x1
y1
x2
y2
x3
y3
ex 3 -------------------------------------------------- */
x1 = 0.0;
= 0.0;
= sqrt(3.0);
= 1.0;
= 1.0;
= sqrt(3.0); */
ex 4 -------------------------------------------------- */
= 0.0;
= 0.0;
= 5.0;
= 2.0;
= 3.0;
= 5.0;
printf("x1 y1 %lg %lg\n", x1, y1);
printf("x2 y2 %lg %lg\n", x2, y2);
printf("x3 y3 %lg %lg\n", x3, y3);
/* a = ((x2-x1)*(y3-y1)-(x3-x1)*(y2-y1)) / 2.0;*/
a = (x2*y3 + x1*y2 + x3*y1 - x2*y1 - x3*y2 - x1*y3) / 2.0;
printf("area %lg\n", a);
}
• これを元に、与えられたメッシュの面積を求めるプログラム例を tri_area.f,
tri_area.c とする。具体的には、ファイル読み込みと、各要素の面積の和
をとるループが追加されている。
• 半径を r = 2 として分割数 m = 3 のとき面積は 12.423 になる。もしこれか
ら大きくずれていたらバグである。デバッグするときは、
1. 入力データがファイルから正しく読み込まれていることを確認する
2. 各要素の面積を書き出して、すべて正の値になっていることを確認する。
3. もし、面積が負になっているものがあったら、先に確認した公式が正し
く導入されていればメッシュデータから三角形を構成するところにバグ
がある。
4. 特に、fortran と C でデータファイルを共有している都合から、データ
ファイルの節点番号は 1 から始まっているが、C では節点座標読み込
み時に節点番号を 0 からつけることになるので、プログラム内の節点
番号と lnods3 に記憶している節点番号は 1 異なっている。
• メッシュ生成プログラムで分割数を変えて数パターン試してみると、分割
数が増えるにしたがって、厳密解に収束していく様子がわかる。例えば、半
径 r = 2 とすると厳密解は 12.566 に対して、分割数 m = 3 のとき面積は
12.423(誤差 1.14%), m = 10 のとき、12.553(誤差 0.103%) となる。
31
!########################################################################
!
tri_area.f
!########################################################################
implicit none
real(8) :: coords(2,10000), pi, area, r, a,
$
x1, y1, x2, y2, x3, y3
integer :: nelem3, nnode, lnods3(3,10000), ielem, i, idum
pi = 3.14159265358979323846264338327950288
area = 0.D0
r = 2.D0
!
!
!
!
!
!
open(10,file=’tri.dat’)
read(10,*) nnode, nelem3
do i = 1, nnode
read(10,*) idum, coords(1:2,i)
write(*,*) idum, coords(1:2,i)
end do
do i = 1, nelem3
read(10,*) idum, lnods3(1:3,i)
write(*,*) idum, lnods3(1:3,i)
end do
close(10)
do ielem = 1, nelem3
x1 = coords(1,lnods3(1,ielem))
y1 = coords(2,lnods3(1,ielem))
x2 = coords(1,lnods3(2,ielem))
y2 = coords(2,lnods3(2,ielem))
x3 = coords(1,lnods3(3,ielem))
y3 = coords(2,lnods3(3,ielem))
write(*,*) x1, y1
write(*,*) x2, y2
write(*,*) x3, y3
a = ((x2-x1)*(y3-y1)-(x3-x1)*(y2-y1)) / 2.D0
a = (x2*y3 + x1*y2 + x3*y1 - x2*y1 - x3*y2 - x1*y3) / 2.D0
write(*,*) a
area = area + a
end do
write(*,*) area*4.D0, pi * r**2
stop
end
/**********************************************************************
tri_area.c
**********************************************************************/
#include <stdio.h>
main()
{
double coords[2][10000], pi, area, r, a, x1, y1, x2, y2, x3, y3;
int nelem3, nnode, lnods3[3][10000], ielem, i, idum;
FILE *fp;
pi = 3.14159265358979323846264338327950288;
area = 0.0;
r = 2.0;
if ((fp = fopen("tri.dat","r")) == NULL)
exit(1);
32
fscanf(fp, "%d %d", &nnode, &nelem3);
for (i=0; i<nnode; i++){
fscanf(fp, "%d %lg %lg", &idum, &coords[0][i], &coords[1][i] );
/*
printf("%d %lg %lg\n", idum, coords[0][i], coords[1][i] );*/
}
for (i=0; i<nelem3; i++){
fscanf(fp, "%d %d %d %d", &idum,
&lnods3[0][i],&lnods3[1][i],&lnods3[2][i]);
/*
printf("%d %d %d %d\n", idum,
lnods3[0][i],lnods3[1][i],lnods3[2][i]);*/
}
for(ielem=0; ielem<nelem3; ielem++){
x1 = coords[0][lnods3[0][ielem]-1];
y1 = coords[1][lnods3[0][ielem]-1];
x2 = coords[0][lnods3[1][ielem]-1];
y2 = coords[1][lnods3[1][ielem]-1];
x3 = coords[0][lnods3[2][ielem]-1];
y3 = coords[1][lnods3[2][ielem]-1];
/*
printf("%lg %lg\n", x1, y1);
printf("%lg %lg\n", x2, y2);
printf("%lg %lg\n", x3, y3);*/
a = ((x2-x1)*(y3-y1)-(x3-x1)*(y2-y1)) / 2.0;
a = (x2*y3 + x1*y2 + x3*y1 - x2*y1 - x3*y2 - x1*y3) / 2.0;
printf("%lg\n", a);
area += a;
}
printf("%lg %lg\n", area*4.0, pi * r*r);
}
33
四角形要素による円の面積
• 次に四角形要素の場合について考える。今回は領域を分割したあとに、さら
に領域内で変数変換を行う。
1 1
A=
dX1 dX2 =
dX1 dX2 =
det[J]dr1 dr2 (19)
V
e
Ve
ここで
e
[J] =
∂X1
∂r1
∂X1
∂r2
∂X2
∂r1
∂X2
∂r2
−1
−1
(20)
さらにこれを Gauss 積分によって計算する。
A=
e
1
−1
1
−1
det[J]dr1 dr2 ≈
2 2
e
Wa Wb det[J(ra , rb )]dr1 dr2 (21)
a=1 b=1
式からわかるように
が 3 つ重なっている。これはプログラムでは 3 重
ループになることに相当する。具体的には、三角形の場合と同じように、最
も外側のループは要素数のループであるが、その内側に、Gauss 積分のルー
プがある。
34
Gauss 積分の実装
• Gauss 積分の実装方法はいろいろあると思うが、ここでは単純にサンプリン
グ点数 nint で示し、配列を用意してサンプリング点の位置 gauss と対応す
る重み weight を記憶する。以下のようにプログラムに直接記入してそのつ
ど関係ないところをコメントアウトすればよい。
1 点積分の場合、
nint = 1 gauss(1) = 0.0
weight(1) = 2.0
2 点積分の場合、
nint = 2 gauss(1) = -sqrt(1.0/3.0)
gauss(2) = sqrt(1.0/3.0)
weight(1) = 1.0
weight(2) = 1.0
3 点積分の場合、
nint = 3 gauss(1) = -sqrt(3.0/5.0)
gauss(2) = 0.0
gauss(3) = sqrt(3.0/5.0)
weight(1) = 5.0/9.0
weight(2) = 8.0/9.0
weight(3) = 5.0/9.0
• 実際にいくつかの関数を積分して比較してみる。
35
1
3 2
=
3x + 2dx = x + 2x
=4
2
−1
−1
1
1
1 3 3 2
1 2
13
x − 3x + 2dx = x − x + 2x
= 4.333 · · ·
=
=
6
2
3
−1 2
−1
1
1
=
2x3 − 4x2 − x − 2dx
5
−1
1
2 4 4 3
1 2
20
= x − x − x − 2x
= − = −6.666 · · ·
4
3
10
3
−1
1
=
−3x4 − −x3 + 2x2 + 5x + 1dx
I1
I2
I3
I4
1
(22)
(23)
(24)
−1
1
3 5 1 4 2 3 5 2
32
= − x − x + x + x +x
= − = 2.1333 · · ·
5
4
3
2
15
−1
π/2
π/2
1
2
2
π
4
I5 =
cos θdθ =
sin θ
cos xdx =
=
2
π
π
−1
−π/2 π
−π/2
(25)
(26)
解析解と、数値積分値を比較すると以下の表のようになる。 2n − 1 次まで
ならば、両者は一致することがわかる。また、多項式ではない cos も、それ
なりによい近似になっていることが確認できる。実際には、解が sin, cos に
なることがわかっているのであれば、その半周期分を 1 要素で表すなどとい
うことはせず、十分細かく要素分割する。
I1
I2
I3
I4
I5
解析解
4
4.333· · ·
-6.666· · ·
2.1333· · ·
1.2732· · ·
1point
4
4
-4
2
2
2point
4
4.333· · ·
-6.666· · ·
2.666· · ·
1.2323· · ·
3point
4
4.333· · ·
-6.666· · ·
2.1333· · ·
1.2741· · ·
• このサンプルコードを gauss.f, gauss.c とする。
!########################################################################
!
gauss.f
!########################################################################
implicit none
real(8) :: gauss(3), weight(3), r, w, f, pi, b
integer :: iint, nint
pi = 3.14159265358979323846264338327950288
!
!
!
nint = 1
gauss(1) = 0.D0
weight(1) = 2.D0
36
!
!
!
!
!
nint = 2
gauss(1) = -dsqrt(1.D0/3.D0)
gauss(2) = dsqrt(1.D0/3.D0)
weight(1) = 1.D0
weight(2) = 1.D0
nint = 3
gauss(1) = -dsqrt(3.D0/5.D0)
gauss(2) = 0.D0
gauss(3) = dsqrt(3.D0/5.D0)
weight(1) = 5.D0/9.D0
weight(2) = 8.D0/9.D0
weight(3) = 5.D0/9.D0
b = 0.D0
do iint = 1, nint
r = gauss(iint)
w = weight(iint)
!
!
!
!
f = 3.D0 * r
+ 2.D0
f = 0.5D0* r**2 - 3.D0 * r
+ 2.D0
f = 2.D0 * r**3 - 4.D0 * r**2 - 0.2D0* r
- 2.D0
f = -3.D0 * r**4 - 1.D0 * r**3 + 2.D0 * r**2 + 5.D0 * r + 1.D0
f = dcos(pi/2.D0 * r)
b = b + w * f
end do
write(*,*) b, 4.D0/pi
stop
end
/**********************************************************************
gauss.c
**********************************************************************/
#include <stdio.h>
main()
{
double gauss[3], weight[3], r, w, f, pi, b;
int iint, nint;
pi = 3.14159265358979323846264338327950288;
nint = 1;
/* gauss[0] = 0.0;
weight[0] = 2.0;*/
/* nint = 2;
gauss[0] = -sqrt(1.0/3.0);
gauss[1] = sqrt(1.0/3.0);
weight[0] = 1.0;
weight[1] = 1.0;*/
nint = 3;
gauss[0] = -sqrt(3.0/5.0);
gauss[1] = 0.0;
gauss[2] = sqrt(3.0/5.0);
weight[0] = 5.0/9.0;
weight[1] = 8.0/9.0;
weight[2] = 5.0/9.0;
37
b = 0.0;
for(iint=0; iint<nint; iint++){
r = gauss[iint];
w = weight[iint];
/*
/*f
/*f
/*f
f =
f = 3.0 * r
+ 2.0;*/
= 0.50* r*r - 3.0 * r
+ 2.0;*/
= 2.0 * r*r*r - 4.0 * r*r - 0.20 * r
- 2.0;*/
= -3.0 * r*r*r*r - 1.0 * r*r*r + 2.0 * r*r + 5.0 * r + 1.0;*/
cos(pi/2.0 * r) ;
b += w * f;
}
printf("%lg %lg", b, 4.0/pi);
}
38
四角形要素による円の面積
•
(n)
∂Xi
∂rj
を求めるための ∂N
は、下式をそのままコーディングするのが結果的
∂rj
にバグも少なく計算機での実行速度も速いようである。+, − が規則的に並ん
でいることに注意すれば確認は比較的容易である。
∂N (1)
∂r1
∂N (2)
∂r1
∂N (3)
∂r1
∂N (4)
∂r1
•
∂Xi
∂rj
1
= − (1 − r2 ),
4
1
(1 − r2 ),
=
4
1
(1 + r2 ),
=
4
1
= − (1 + r2 ),
4
∂N (1)
∂r2
∂N (2)
∂r2
∂N (3)
∂r2
∂N (3)
∂r2
1
= − (1 − r1 )
4
1
= − (1 + r1 )
4
1
(1 + r1 )
=
4
1
=
(1 − r1 )
4
(27)
(28)
(29)
(30)
もストレートにコーディングするほうが簡単である。
∂ (1) (1)
∂Xi
(2)
(3)
(4)
=
N Xi + N (2) Xi + N (3) Xi + N (4) Xi
∂rj
∂rj
∂N (1) (1) ∂N (2) (2) ∂N (3) (3) ∂N (4) (4)
Xi +
Xi +
Xj +
Xj
=
∂rj
∂rj
∂ri
∂ri
(31)
(32)
• 実際のプログラム例を quad_area.f, quad_area.c とする。面積を求める場
合は、四角形要素を用いても、三角形要素と同じ計算結果になることを確認
できる。また、たまたま今回の問題では積分点数によらず同じ解になる。例
えば、半径 r = 2 とすると厳密解は 12.566 に対して、分割数 m = 3 のとき
面積は 12.423(誤差 1.14%), m = 10 のとき、12.553(誤差 0.103%) となる。
• もしこの値から大幅に外れるときは、バグである。デバッグの注意点は三角
形要素と同じであるが、
1. 各要素の面積が正、あるいは、合理的な値になっているか、すなわち、
桁違いに大きかったり小さかったりしないことを確認する。
2. これが正しければ、各要素の面積を加算するところなど単純なところに
バグがある可能性が高い。
3. 面積がおかしいときは、まず、Gauus 積分の 2 重ループの中で r1, r1 , w1 , w2
の値が配列 gauss, weight から正しく代入されているか確認する。
4. Gauss 積分単独なら正しく動作する場合、入力データを (−1, −1), (1, −1), (1, 1), (−1, 1)
の正方形 1 要素をデータとして与えると、正しくコーディングされていれ
ば、積分点によらず、ヤコビマトリックスは単位行列でかつ、det[J] = 1
であることを利用して、デバッグできる。
39
5. また、三角形同様に、fortran と C でデータファイルを共有している都
合から、データファイルの節点番号は 1 から始まっていることに対応
しているか確認する。C では節点座標読み込み時に節点番号を 0 から
つけることになるので、プログラム内の節点番号と lnods に記憶して
いる節点番号は 1 異なっている。
• メッシュ生成プログラムで分割数を変えて数パターン試してみると、分割
数が増えるにしたがって、厳密解に収束していく様子がわかる。例えば、半
径 r = 2 とすると厳密解は 12.566 に対して、分割数 m = 3 のとき面積は
12.423(誤差 1.14%), m = 10 のとき、12.553(誤差 0.103%) となる。
!########################################################################
!
quad_area.f
!########################################################################
implicit none
real(8) :: coords(2,10000), pi, area, r, gauss(3), weight(3),
$
r1, r2, w1, w2, jacobi(2,2), detj, n(4), dndr(2,4)
integer :: nelem, nnode, lnods(4,10000), ielem, i, idum,
$
iint1, iint2, nint
!
!
!
nint = 1
gauss(1) = 0.D0
weight(1) = 2.D0
!
!
!
!
!
nint = 2
gauss(1) = -dsqrt(1.D0/3.D0)
gauss(2) = dsqrt(1.D0/3.D0)
weight(1) = 1.D0
weight(2) = 1.D0
nint = 3
gauss(1) = -dsqrt(3.D0/5.D0)
gauss(2) = 0.D0
gauss(3) = dsqrt(3.D0/5.D0)
weight(1) = 5.D0/9.D0
weight(2) = 8.D0/9.D0
weight(3) = 5.D0/9.D0
pi = 3.14159265358979323846264338327950288
area = 0.D0
r = 2.D0
!
!
open(10,file=’quad.dat’)
read(10,*) nnode, nelem
do i = 1, nnode
read(10,*) idum, coords(1:2,i)
write(*,*) idum, coords(1:2,i)
end do
do i = 1, nelem
read(10,*) idum, lnods(1:4,i)
write(*,*) idum, lnods(1:4,i)
end do
close(10)
do ielem = 1, nelem
do iint1 = 1, nint
do iint2 = 1, nint
r1 = gauss(iint1)
r2 = gauss(iint2)
40
w1 = weight(iint1)
w2 = weight(iint2)
dndr(1,1)
dndr(1,2)
dndr(1,3)
dndr(1,4)
= -1.D0/4.D0 * (1-r2)
= 1.D0/4.D0 * (1-r2)
= 1.D0/4.D0 * (1+r2)
= -1.D0/4.D0 * (1+r2)
dndr(2,1)
dndr(2,2)
dndr(2,3)
dndr(2,4)
= -1.D0/4.D0 * (1-r1)
= -1.D0/4.D0 * (1+r1)
= 1.D0/4.D0 * (1+r1)
= 1.D0/4.D0 * (1-r1)
$
jacobi(1:2,1:2)=0.D0
do i = 1, 4
jacobi(1,1) = jacobi(1,1)
+ dndr(1,i)*coords(1,lnods(i,ielem))
$
jacobi(1,2) = jacobi(1,2)
+ dndr(1,i)*coords(2,lnods(i,ielem))
$
jacobi(2,1) = jacobi(2,1)
+ dndr(2,i)*coords(1,lnods(i,ielem))
$
jacobi(2,2) = jacobi(2,2)
+ dndr(2,i)*coords(2,lnods(i,ielem))
end do
$
detj = jacobi(1,1)*jacobi(2,2)
- jacobi(1,2)*jacobi(2,1)
area = area + w1 * w2 * detj
end do
end do
end do
write(*,*) area*4.D0, pi * r**2
stop
end
/**********************************************************************
quad_area.c
**********************************************************************/
#include <stdio.h>
main()
{
double coords[2][10000], pi, area, r, gauss[3], weight[3],
r1, r2, w1, w2, jacobi[2][2], detj, n[4], dndr[2][4];
int nelem, nnode, lnods[4][10000], ielem, i,j, idum, iint1, iint2, nint;
FILE *fp;
/*nint = 1;
gauss[0] = 0.0;
weight[0] = 2.0;*/
/* nint = 2;
gauss[0] = -sqrt(1.0/3.0);
gauss[1] = sqrt(1.0/3.0);
weight[0] = 1.0;
weight[1] = 1.0;*/
nint = 3;
gauss[0] = -sqrt(3.0/5.0);
41
gauss[1] =
gauss[2] =
weight[0] =
weight[1] =
weight[2] =
0.0;
sqrt(3.0/5.0);
5.0/9.0;
8.0/9.0;
5.0/9.0;
pi = 3.14159265358979323846264338327950288;
area = 0.0;
r = 2.0;
if ((fp = fopen("quad.dat","r")) == NULL)
exit(1);
fscanf(fp, "%d %d", &nnode, &nelem);
for (i=0; i<nnode; i++){
fscanf(fp, "%d %lg %lg", &idum, &coords[0][i], &coords[1][i] );
/*
printf("%d %lg %lg\n", idum, coords[0][i], coords[1][i] );*/
}
for (i=0; i<nelem; i++){
fscanf(fp, "%d %d %d %d %d", &idum,
&lnods[0][i],&lnods[1][i],&lnods[2][i],&lnods[3][i]);
/*
printf("%d %d %d %d %d\n", idum,
lnods[0][i],lnods[1][i],lnods[2][i],lnods[3][i]);*/
}
for(ielem=0; ielem<nelem; ielem++){
for(iint1=0; iint1<nint; iint1++){
for(iint2=0; iint2<nint; iint2++){
r1 = gauss[iint1];
r2 = gauss[iint2];
w1 = weight[iint1];
w2 = weight[iint2];
dndr[0][0]
dndr[0][1]
dndr[0][2]
dndr[0][3]
= -1.0/4.0 * (1.0-r2);
= 1.0/4.0 * (1.0-r2);
= 1.0/4.0 * (1.0+r2);
= -1.0/4.0 * (1.0+r2);
dndr[1][0]
dndr[1][1]
dndr[1][2]
dndr[1][3]
= -1.0/4.0 * (1.0-r1);
= -1.0/4.0 * (1.0+r1);
= 1.0/4.0 * (1.0+r1);
= 1.0/4.0 * (1.0-r1);
for(i=0; i<2; i++){
for(j=0; j<2; j++){
jacobi[i][j] = 0.0;
}
}
for (i=0; i<4; i++){
jacobi[0][0] += dndr[0][i] * coords[0][lnods[i][ielem]-1];
jacobi[0][1] += dndr[0][i] * coords[1][lnods[i][ielem]-1];
jacobi[1][0] += dndr[1][i] * coords[0][lnods[i][ielem]-1];
jacobi[1][1] += dndr[1][i] * coords[1][lnods[i][ielem]-1];
}
detj = jacobi[0][0] * jacobi[1][1]
- jacobi[0][1] * jacobi[1][0];
area = area + w1 * w2 * detj;
}
}
}
printf("%lg, %lg\n", area*4.0, pi*r*r);
}
42
三角形要素による球の体積
• 円の面積を求めるプログラムをベースにして球の体積を求めるプログラムを
作成する。実は有限要素法のプログラムは円の面積のプログラムをベースに
して作成できるので、球の体積は本筋とは関係ないが、要素種類と積分との
関連性のイメージを持ってもらうためにあえて取り上げる。
• 領域積分の形式で、球の体積 W を求める式は以下のとおりである。
W =
r 2 − X12 − X22 dX1 dX2
(33)
V
• 三角形要素を使った有限要素法の手順で球の体積を求める場合、三角形領域
の内部では被積分関数は定数になることを考慮して、積分は定数面積の積に
なるということをそのまま用いることにする。すなわち、三角形の重心を G
として
h=
r 2 − G21 − G22
(34)
により、三角柱の高さを求め、個々の体積の和を取るというプロセスになる。
すなわち三角形の面積を Ae として
2
2
2
W =
r − X1 − X2 dX1 dX2 =
r 2 − X12 − X22 dX1 dX2
V
=
e
Ve
(35)
e
A r 2 − G21 − G22
(36)
e
• このサンプルプログラムを tri_int.f, tri_int.c とする。ただし、体積 W
は変数 vol, 重心の座標は g1, g2, 重心における三角柱の高さは h として
いる。
• 半径 r = 2 とすると、分割数 m = 3 のとき体積は 33.87 になる。もしこの
値から大きくずれていたらバグがある。データファイルから読み込んだ円の
メッシュの面積が正しく求められているのであれば、ほとんど流用できるの
で、バグの可能性は重心の計算、三角柱の高さの計算などに限られている。
• 面積同様メッシュ生成プログラムで分割数を変えて数パターン試してみると、
分割数が増えるにしたがって、厳密解に収束していく様子がわかる。例えば、
半径 r = 2 とすると厳密解は 33.51 に対して、分割数 m = 3 のとき面積は
33.87(誤差 0.865%), m = 10 のとき、33.58(誤差 0.209%) となる。
43
!########################################################################
!
tri_int.f
!########################################################################
implicit none
real(8) :: coords(2,10000), pi, vol, r, a, h, g1, g2,
$
x1, y1, x2, y2, x3, y3
integer :: nelem, nnode, lnods(3,10000), ielem, i, idum
pi = 3.14159265358979323846264338327950288
vol = 0.D0
r = 2.D0
!
!
open(10,file=’tri.dat’)
read(10,*) nnode, nelem
do i = 1, nnode
read(10,*) idum, coords(1:2,i)
write(*,*) idum, coords(1:2,i)
end do
do i = 1, nelem
read(10,*) idum, lnods(1:3,i)
write(*,*) idum, lnods(1:3,i)
end do
close(10)
do ielem = 1, nelem
x1 = coords(1,lnods(1,ielem))
y1 = coords(2,lnods(1,ielem))
x2 = coords(1,lnods(2,ielem))
y2 = coords(2,lnods(2,ielem))
x3 = coords(1,lnods(3,ielem))
y3 = coords(2,lnods(3,ielem))
write(*,*) x1, y1
write(*,*) x2, y2
write(*,*) x3, y3
g1 = (x1 + x2 + x3) / 3.D0
g2 = (y1 + y2 + y3) / 3.D0
write(*,*) ielem, g1, g2
h = dsqrt(r**2 - g1**2 - g2**2)
a = (x2*y3 + x1*y2 + x3*y1 - x2*y1 - x3*y2 - x1*y3) / 2.D0
write(*,*) h, a
vol = vol + h * a
end do
write(*,*) vol*8.D0, 4.D0 * pi * r**3 / 3
stop
end
/**********************************************************************
tri_int.c
**********************************************************************/
#include <stdio.h>
main()
{
double coords[2][10000], pi, vol, r, a, h, g1, g2, x1, y1, x2, y2, x3, y3;
int nelem3, nnode, lnods3[3][10000], ielem, i, idum;
FILE *fp;
pi = 3.14159265358979323846264338327950288;
vol = 0.0;
44
r = 2.0;
if ((fp = fopen("tri.dat","r")) == NULL)
exit(1);
fscanf(fp, "%d %d", &nnode, &nelem3);
for (i=0; i<nnode; i++){
fscanf(fp, "%d %lg %lg", &idum, &coords[0][i], &coords[1][i] );
/*
printf("%d %lg %lg\n", idum, coords[0][i], coords[1][i] );*/
}
for (i=0; i<nelem3; i++){
fscanf(fp, "%d %d %d %d", &idum,
&lnods3[0][i],&lnods3[1][i],&lnods3[2][i]);
/*
printf("%d %d %d %d\n", idum,
lnods3[0][i],lnods3[1][i],lnods3[2][i]);*/
}
for(ielem=0; ielem<nelem3; ielem++){
x1 = coords[0][lnods3[0][ielem]-1];
y1 = coords[1][lnods3[0][ielem]-1];
x2 = coords[0][lnods3[1][ielem]-1];
y2 = coords[1][lnods3[1][ielem]-1];
x3 = coords[0][lnods3[2][ielem]-1];
y3 = coords[1][lnods3[2][ielem]-1];
/*
printf("%lg %lg\n", x1, y1);
printf("%lg %lg\n", x2, y2);
printf("%lg %lg\n", x3, y3);*/
g1 = (x1 + x2 + x3) / 3.0;
g2 = (y1 + y2 + y3) / 3.0;
h = sqrt(r*r - g1*g1 - g2*g2);
a = (x2*y3 + x1*y2 + x3*y1 - x2*y1 - x3*y2 - x1*y3) / 2.0;
/*
printf("%lg\n", a);*/
vol += h * a;
}
printf("%lg %lg\n", vol*8.0, 4.0*pi * r*r*r/3.0);
}
45
四角形要素による球の体積
• 四角形要素を用いる場合、領域積分の形式で、球の体積 W を求める式を変
数変換し、さらに Gauss 積分を適応する。
W =
r 2 − X12 − X22 dX1 dX2
(37)
V
=
r 2 − X12 − X22 dX1 dX2
(38)
e
=
e
≈
Ve
1
−1
1
−1
2 2
e
r 2 − X12 − X22 det[J]dr1 dr2
Wa Wb
(39)
r 2 − X12 − X22 det[J(ra , rb )]dr1 dr2
(40)
a=1 b=1
面積を求めるプログラムは被積分関数は det[J] だけであったが、今度は、
r 2 − X12 − X22 が加わる。サンプリング点における X1 , X2 の値は、
(1)
(2)
(3)
(4)
(1)
(2)
(3)
(4)
X1 = N (1) X1 + N (2) X1 + N (3) X1 + N (4) X1
X2 = N (1) X2 + N (2) X2 + N (3) X2 + N (4) X2
(41)
(42)
(43)
によって求められる。なお、補間関数は以下のとおりである。
1
N (1) = (1 − r1 )(1 − r2 )
4
1
N (2) = (1 + r1 )(1 − r2 )
4
1
N (3) = (1 + r1 )(1 + r2 )
4
1
N (4) = (1 − r1 )(1 + r2 )
4
(44)
(45)
(46)
(47)
• これを元に作成したサンプルプログラムを quad_int.f, quad_int.c とする。
• Gauss 積分のサンプリング点数 3、半径 r = 2 とすると、分割数 m = 3 の
とき体積は 33.48 になる。もしこの値から大きくずれていたらバグがある。
データファイルから読み込んだ円のメッシュの面積が正しく求められている
のであれば、三角形要素の場合と同様に、プログラムは面積を求めるプログ
ラムへの追加でコーディングできるので、バグの可能性は新たにコーディン
グした、補間関数 N 、座標 X1 , X2 の計算、高さを求めるところなどに限定
される。
46
• 面積同様メッシュ生成プログラムで分割数を変えて数パターン試してみると、
分割数が増えるにしたがって、厳密解に収束していく様子がわかる。また、
例えば、半径 r = 2 とすると厳密解は 33.51 に対して、分割数 m = 3 のとき
面積は 33.48(誤差 0.865%), m = 10 のとき、33.58(誤差 0.209%) となる。
• 実際の有限要素法のプログラムも、最初からすべての定式化をコーディング
すると確実に破綻する。今回のように「きりわけ」ができるところを探して、
ひとつずつ積み上げていくと、中間のプログラムをたくさん作ることになる
ので、コーディングとしては多少遠回りをしている感じがするものの、結果
的にはデバッグの時間を大幅に短縮できることになる。
!########################################################################
!
quad_int.f
!########################################################################
implicit none
real(8) :: coords(2,10000), pi, vol, r, f,
$
gauss(3), weight(3),
$
r1, r2, w1, w2, jacobi(2,2), detj,x,y,
$
n(4), dndr(2,4), s
integer :: nelem, nnode, lnods(4,10000), ielem, i, idum,
$
iint1, iint2, nint
!
!
!
nint = 1
gauss(1) = 0.D0
weight(1) = 2.D0
nint = 2
gauss(1) = -dsqrt(1.D0/3.D0)
gauss(2) = dsqrt(1.D0/3.D0)
weight(1) = 1.D0
weight(2) = 1.D0
nint = 3
gauss(1) = -dsqrt(3.D0/5.D0)
gauss(2) = 0.D0
gauss(3) = dsqrt(3.D0/5.D0)
weight(1) = 5.D0/9.D0
weight(2) = 8.D0/9.D0
weight(3) = 5.D0/9.D0
pi = 3.14159265358979323846264338327950288
vol = 0.D0
s = 0.D0
r = 2.D0
!
!
open(10,file=’quad.dat’)
read(10,*) nnode, nelem
do i = 1, nnode
read(10,*) idum, coords(1:2,i)
write(*,*) idum, coords(1:2,i)
end do
do i = 1, nelem
read(10,*) idum, lnods(1:4,i)
write(*,*) idum, lnods(1:4,i)
end do
close(10)
do ielem = 1, nelem
do iint1 = 1, nint
do iint2 = 1, nint
47
r1
r2
w1
w2
=
=
=
=
gauss(iint1)
gauss(iint2)
weight(iint1)
weight(iint2)
dndr(1,1)
dndr(1,2)
dndr(1,3)
dndr(1,4)
= -1.D0/4.D0 * (1-r2)
= 1.D0/4.D0 * (1-r2)
= 1.D0/4.D0 * (1+r2)
= -1.D0/4.D0 * (1+r2)
dndr(2,1)
dndr(2,2)
dndr(2,3)
dndr(2,4)
= -1.D0/4.D0 * (1-r1)
= -1.D0/4.D0 * (1+r1)
= 1.D0/4.D0 * (1+r1)
= 1.D0/4.D0 * (1-r1)
$
jacobi(1:2,1:2) = 0.D0
do i = 1, 4
jacobi(1,1) = jacobi(1,1)
+ dndr(1,i)*coords(1,lnods(i,ielem))
$
jacobi(1,2) = jacobi(1,2)
+ dndr(1,i)*coords(2,lnods(i,ielem))
$
jacobi(2,1) = jacobi(2,1)
+ dndr(2,i)*coords(1,lnods(i,ielem))
$
jacobi(2,2) = jacobi(2,2)
+ dndr(2,i)*coords(2,lnods(i,ielem))
end do
$
detj = jacobi(1,1)*jacobi(2,2)
- jacobi(1,2)*jacobi(2,1)
n(1)
n(2)
n(3)
n(4)
=
=
=
=
1.D0/4.D0
1.D0/4.D0
1.D0/4.D0
1.D0/4.D0
*
*
*
*
(1.D0
(1.D0
(1.D0
(1.D0
+
+
-
r1)
r1)
r1)
r1)
*
*
*
*
(1.D0
(1.D0
(1.D0
(1.D0
$
$
$
x = n(1) * coords(1,lnods(1,ielem))
+ n(2) * coords(1,lnods(2,ielem))
+ n(3) * coords(1,lnods(3,ielem))
+ n(4) * coords(1,lnods(4,ielem))
$
$
$
y = n(1) * coords(2,lnods(1,ielem))
+ n(2) * coords(2,lnods(2,ielem))
+ n(3) * coords(2,lnods(3,ielem))
+ n(4) * coords(2,lnods(4,ielem))
+
+
r2)
r2)
r2)
r2)
f = dsqrt(r**2 - x**2 - y**2)
vol = vol + f * w1 * w2 * detj
end do
end do
end do
write(*,*) vol*8.D0, 4.D0 * pi * r**3 / 3
stop
end
/**********************************************************************
quad_int.c
**********************************************************************/
#include <stdio.h>
48
main()
{
double coords[2][10000], pi, vol, r, f, gauss[3], weight[3],
r1, r2, w1, w2, jacobi[2][2], detj, x, y, n[4], dndr[2][4];
int nelem, nnode, lnods[4][10000], ielem, i,j, idum, iint1, iint2, nint;
FILE *fp;
/*nint = 1;
gauss[0] = 0.0;
weight[0] = 2.0;*/
/* nint = 2;
gauss[0] = -sqrt(1.0/3.0);
gauss[1] = sqrt(1.0/3.0);
weight[0] = 1.0;
weight[1] = 1.0;*/
nint = 3;
gauss[0] = -sqrt(3.0/5.0);
gauss[1] = 0.0;
gauss[2] = sqrt(3.0/5.0);
weight[0] = 5.0/9.0;
weight[1] = 8.0/9.0;
weight[2] = 5.0/9.0;
pi = 3.14159265358979323846264338327950288;
vol = 0.0;
r = 2.0;
if ((fp = fopen("quad.dat","r")) == NULL)
exit(1);
fscanf(fp, "%d %d", &nnode, &nelem);
for (i=0; i<nnode; i++){
fscanf(fp, "%d %lg %lg", &idum, &coords[0][i], &coords[1][i] );
/*
printf("%d %lg %lg\n", idum, coords[0][i], coords[1][i] );*/
}
for (i=0; i<nelem; i++){
fscanf(fp, "%d %d %d %d %d", &idum,
&lnods[0][i],&lnods[1][i],&lnods[2][i],&lnods[3][i]);
/*
printf("%d %d %d %d %d\n", idum,
lnods[0][i],lnods[1][i],lnods[2][i],lnods[3][i]);*/
}
for(ielem=0; ielem<nelem; ielem++){
for(iint1=0; iint1<nint; iint1++){
for(iint2=0; iint2<nint; iint2++){
r1 = gauss[iint1];
r2 = gauss[iint2];
w1 = weight[iint1];
w2 = weight[iint2];
dndr[0][0]
dndr[0][1]
dndr[0][2]
dndr[0][3]
= -1.0/4.0 * (1.0-r2);
= 1.0/4.0 * (1.0-r2);
= 1.0/4.0 * (1.0+r2);
= -1.0/4.0 * (1.0+r2);
dndr[1][0]
dndr[1][1]
dndr[1][2]
dndr[1][3]
= -1.0/4.0 * (1.0-r1);
= -1.0/4.0 * (1.0+r1);
= 1.0/4.0 * (1.0+r1);
= 1.0/4.0 * (1.0-r1);
for(i=0; i<2; i++){
for(j=0; j<2; j++){
jacobi[i][j] = 0.0;
}
49
}
for (i=0; i<4; i++){
jacobi[0][0] += dndr[0][i] * coords[0][lnods[i][ielem]-1];
jacobi[0][1] += dndr[0][i] * coords[1][lnods[i][ielem]-1];
jacobi[1][0] += dndr[1][i] * coords[0][lnods[i][ielem]-1];
jacobi[1][1] += dndr[1][i] * coords[1][lnods[i][ielem]-1];
}
detj = jacobi[0][0] * jacobi[1][1]
- jacobi[0][1] * jacobi[1][0];
n[0]
n[1]
n[2]
n[3]
=
=
=
=
1.0/4.0
1.0/4.0
1.0/4.0
1.0/4.0
*
*
*
*
(1.0
(1.0
(1.0
(1.0
+
+
-
r1)
r1)
r1)
r1)
*
*
*
*
(1.0
(1.0
(1.0
(1.0
+
+
r2);
r2);
r2);
r2);
x = n[0] * coords[0][lnods[0][ielem]-1]
+ n[1] * coords[0][lnods[1][ielem]-1]
+ n[2] * coords[0][lnods[2][ielem]-1]
+ n[3] * coords[0][lnods[3][ielem]-1];
y = n[0] * coords[1][lnods[0][ielem]-1]
+ n[1] * coords[1][lnods[1][ielem]-1]
+ n[2] * coords[1][lnods[2][ielem]-1]
+ n[3] * coords[1][lnods[3][ielem]-1];
f = sqrt(r*r - x*x - y*y);
vol += f * w1 * w2 * detj;
}
}
}
printf("%lg, %lg\n", vol*8.0, 4.0*pi*r*r*r/3.0);
}
50