第7章有限要素法のプログラミング - 名古屋大学

はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
.
.
第 7 章有限要素法のプログラミング
畔上秀幸
名古屋大学情報科学研究科複雑系科学専攻
December 18, 2014
.
.
.
.
.
.
1 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.1 はじめに
(目標) 2 次元 Poisson 問題を有限要素法で解くためのプログラムの
構成について理解する．
.
.
.
.
.
.
2 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.2 有限要素法プログラムの構成
有限要素法プログラムの具体例をみてみよう．授業のウェブページに
あるファイル femfp.c は，テキスト [1] のプログラムを基にして作成さ
れた．図のような関数で構成されている．
main( )
input( )
assem( )
ecm( )
solve( )
gs{solve( )
f( )
output( )
.
.
.
.
.
.
3 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.3 ファイル入出力関数
.
fopen()
.
#include <stdio.h>
FILE *fopen(char *fname, char *mode);
• fname のファイルを mode (下表) にしたがって開く．
• 返却値は，成功したときファイルの番号，失敗したとき NULL と
なる．
mode
r
r+
w
w+
a
a
.
作用
テキストの読み込み
テキストの読み書き
テキストの書き込み
テキストの読み書き
テキストの追加書き込み
テキストの読み書き
mode
rb
rb+
wb
wb+
ab
ab+
作用
バイナリの読み込み
バイナリファイルの読み書き
バイナリの書き込み
バイナリの読み書き
バイナリの追加書き込み
バイナリの読み書き
.
.
.
.
.
.
4 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.3 ファイル入出力関数 (cnt.)
.
fclose()
.
#include <stdio.h>
int fclose(FILE *fp);
• fopen で開いたファイル fp を閉じる．
. • 返却値は，成功したとき 0，失敗したとき EOF となる．
.
fprintf()
.
include <stdio.h>
int fprintf(FILE *fp, char *fmt,...);
• ファイル fp に対して，書式付き文字列 fmt を書き込む．
• 返却値は，成功したとき出力した文字数，失敗したときは負の数と
.
なる．
.
.
.
.
.
.
5 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.3 ファイル入出力関数 (cnt.)
.
fputs()
.
#include <stdio.h>
int fputs(char *string, FILE *fp);
• ファイル fp に対して，文字列 string を書き込む．
. • 返却値は，成功したとき 0，失敗したとき EOF となる．
.
fgets()
.
#include <stdio.h>
char *fgets(char *string, int n, FILE *fp);
• ファイル fp に対して，文字列 string を読み込む．
. • 返却値は，成功したとき 0，失敗したとき EOF となる．
.
.
.
.
.
.
6 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.4 有限要素法プログラム
具体的な記述をみてみよう．
.
main()
.
/* The Finite Element Method for the Poisson Equation */
/* Full Matrix Version */
/* The Gauss Elimination Method */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define NDIM1 200
#define
NDIM2 400
.
.
.
.
.
.
.
7 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.4 有限要素法プログラム (cnt.)
.
void input(int *nnode, int *nelmt, int *nbc,
int ielmt[][3], double x[],
double y[], int ibc[]);
void assem(int nnode, int nelmt, int nbc, int ielmt[][3],
double x[], double y[], int ibc[],
double am[][NDIM1], double fm[]);
void ecm(int ielmt[][3], double x[], double y[],
int ie, double ae[][3], double fe[]);
void solve(int m, double a[][NDIM1], double f[]);
void gs_solve(int m, double a[][NDIM1], double f[]);
void output(double fm[], int nnode, double x[], double y[]);
double
f(double x, double y);
.
.
.
.
.
.
.
8 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.4 有限要素法プログラム (cnt.)
.
/* Main Program */
int main()
{
int nnode, nelmt, nbc, ielmt[NDIM2][3], ibc[NDIM1];
double am[NDIM1][NDIM1], fm[NDIM1], x[NDIM1], y[NDIM1];
input(&nnode, &nelmt, &nbc, ielmt, x, y, ibc);
assem(nnode, nelmt, nbc, ielmt, x, y, ibc, am, fm);
solve(nnode, am, fm);
//gs_solve(nnode, am, fm);
output(fm, nnode, x, y);
return(EXIT_SUCCESS);
}
.
テキスト [1] からの修正:
• void output(double fm[], int nnode, double x[], double y[]);
• output(fm, nnode, x, y);
.
.
.
.
.
.
9 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.4 有限要素法プログラム (cnt.)
.
input()
.
/* Input */
void input(int *nnode, int *nelmt, int *nbc,
int ielmt[][3], double x[], double y[], int ibc[])
{
int j;
FILE *fpin, *fpout;
.
if((fpin=fopen("indata.txt","r")) == NULL){
printf("入力データファイルが見つかりません．\n");
exit(EXIT_FAILURE);
}
if((fpout=fopen("outdata.txt","w")) == NULL){
fclose(fpin);
printf("出力データファイルが作成できません．\n");
exit(EXIT_FAILURE);
}
.
.
.
.
.
.
10 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.4 有限要素法プログラム (cnt.)
.
.
/* Input of Data */
printf("Input (nnode, nelmt, nbc): ");
fscanf(fpin, "%d%d%d", nnode, nelmt, nbc);
printf("\nInput(x,y) of node i: \n");
for(j=0; j<*nnode; j++){
printf("
for i=%d: \n", j+1);
fscanf(fpin, "%lg%lg", &(x[j]), &(y[j]));
}
printf("Input(1st, 2nd, 3rd) nodes of element i: \n");
for(j=0; j<*nelmt; j++){
printf("
for i=%d: \n", j+1);
fscanf(fpin, "%d%d%d", &(ielmt[j][0]), &(ielmt[j][1]),
&(ielmt[j][2]));
}
.
.
.
.
.
.
11 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.4 有限要素法プログラム (cnt.)
.
.
if(*nbc>0){
printf("Input nodes with zero Dirichlet data");
printf(" for i=1 to %d: ", *nbc);
for(j=0; j<*nbc; j++) fscanf(fpin, "%d", &(ibc[j]));
}
printf("\n");
printf("Nodes of elements \n");
for(j=0; j<2; j++) printf("
i
i1 i2 i3 ");
printf("\n");
for(j=0; j<*nelmt; j+=2){
printf("*%4d*%4d%4d%4d ", j+1, ielmt[j][0],
ielmt[j][1], ielmt[j][2]);
if(j<*nelmt-1)
printf("*%4d*%4d%4d%4d", j+2, ielmt[j+1][0],
ielmt[j+1][1], ielmt[j+1][2]);
printf("\n");
}
.
.
.
.
.
.
12 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.4 有限要素法プログラム (cnt.)
.
if(*nbc>0){
printf("Nodes with zero Dirichlet data \n");
for(j=0; j<*nbc; j++) printf("%4d", ibc[j]);
printf("\n");
}
}
.
テキスト [1] からの修正:
• FILE *fpin, *fpout; ...
• 入力ファイル名: indata.txt, 出力ファイル名: outdata.txt
• fscanf(fpin, ...
.
.
.
.
.
.
13 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.4 有限要素法プログラム (cnt.)
.
assem()
.
/* The Direct Stiffness Method */
void assem(int nnode, int nelmt, int nbc, int ielmt[][3],
double x[], double y[], int ibc[],
double am[NDIM1][NDIM1], double fm[])
{
int i, j, ie, ii, jj;
double ae[3][3], fe[3];
.
/* Initial Clear */
for(i=0; i<nnode; i++){
fm[i]=0.0;
for(j=0; j<nnode; j++) am[i][j]=0.0;
}
.
.
.
.
.
.
14 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.4 有限要素法プログラム (cnt.)
.
.
/* Assemblage of Total Matrix and Vector */
for(ie=0; ie<nelmt; ie++){
ecm(ielmt, x, y, ie, ae, fe);
for(i=0; i<3; i++){
ii=ielmt[ie][i]-1;
fm[ii]+=fe[i];
for(j=0; j<3; j++){
jj=ielmt[ie][j]-1;
am[ii][jj]+=ae[i][j];
}
}
}
.
.
.
.
.
.
15 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.4 有限要素法プログラム (cnt.)
.
/* The Homogeneous Dirichlet Condition */
if (nbc!=0){
for(i=0;i<nbc; i++){
ii=ibc[i]-1;
fm[ii]=0.0;
for(j=0; j<nnode;j++){
am[ii][j]=0.0;
am[j][ii]=0.0;
}
am[ii][ii]=1.0;
}
}
}
.
.
.
.
.
.
.
16 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.4 有限要素法プログラム (cnt.)
.
ecm()
.
/* Element Coeffcient Matrix and Free Vector */
void ecm(int ielmt[][3], double x[], double y[],
int ie, double ae[][3], double fe[])
{
int i, j, k;
double d, s, xe[3], ye[3], b[3], c[3];
.
for(i=0; i<3;i++){
j=ielmt[ie][i]-1;
xe[i]=x[j];
ye[i]=y[j];
}
d=xe[0]*(ye[1]-ye[2])+xe[1]*(ye[2]-ye[0])
+xe[2]*(ye[0]-ye[1]);
s=fabs(d)/2.0;
.
.
.
.
.
.
17 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.4 有限要素法プログラム (cnt.)
.
/* Calculation of Element Coefficient Matrix */
for(i=0; i<3; i++){
j=i+1; if(i==2) j=0;
k=i-1; if(i==0) k=2;
b[i]=(ye[j]-ye[k])/d;
c[i]=(xe[k]-xe[j])/d;
}
for(i=0; i<3; i++){
for(j=0; j<3; j++) ae[i][j]=s*(b[j]*b[i]+c[j]*c[i]);
}
/* Calculation of Element Free Vector */
for(i=0; i<3; i++) fe[i]=s*f(xe[i], ye[i])/3.0;
}
.
.
.
.
.
.
.
18 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.4 有限要素法プログラム (cnt.)
.
solve()
.
/* The Gauss Elimination Method */
void solve(int m, double a[][NDIM1], double f[])
/* m=number of unknowns */
{
int i, j, k;
double aa;
.
/* Forward Elimination */
for(i=0; i<m-1; i++){
for(j=i+1; j<m; j++){
aa=a[j][i]/a[i][i];
f[j]-=aa*f[i];
for(k=i+1; k<m; k++) a[j][k]-=aa*a[i][k];
}
}
.
.
.
.
.
.
19 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.4 有限要素法プログラム (cnt.)
.
/* Backward Substitution */
f[m-1]/=a[m-1][m-1];
for(i=m-2; i>=0; i--){
for(j=i+1; j<m; j++) f[i]-=a[i][j]*f[j];
f[i]/=a[i][i];
}
}
.
.
.
.
.
.
.
20 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.4 有限要素法プログラム (cnt.)
テキスト [1] からの修正:
gs solve() を追加
.
gs solve()
.
/* The Gauss Seidel Method */
void gs_solve(int m, double a[][NDIM1], double f[])
{
int i, j, k;
double aa = 0.0;
double u[NDIM1],u_old[NDIM1];
double max_delta=0.0, max_delta_old;
double max_u=0.0;
int count = 0;
for(i=0; i<m; i++){
u[i] = 0.0;
}
.
.
.
.
.
.
.
21 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.4 有限要素法プログラム (cnt.)
.
do{
max_delta_old = max_delta;
for(i=0; i<m; i++){
u_old[i] = u[i];
}
.
for(i=0; i<m; i++){
aa = 0.0;
for(j=0; j<m; j++){
if(i != j){
aa += a[i][j]*u[j];
}
}
u[i]=(f[i] - aa)/a[i][i];
}
.
.
.
.
.
.
22 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.4 有限要素法プログラム (cnt.)
.
max_delta = fabs(u[0] - u_old[0]);
for(i=1; i<m; i++){
if(max_delta < fabs(u[i] - u_old[i])){
max_delta = fabs(u[i] - u_old[i]);
}
}
.
max_u = fabs(u[0]);
for(i=1; i<m; i++){
if(max_u < fabs(u[i])){
max_u = fabs(u[i]);
}
}
.
.
.
.
.
.
23 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.4 有限要素法プログラム (cnt.)
.
count++;
}while(fabs((max_delta - max_delta_old)/max_u) >= 1.0E-5);
printf("count = %d\n",count);
for(i=0; i<m; i++){
f[i] = u[i];
}
}
.
.
.
.
.
.
.
24 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.4 有限要素法プログラム (cnt.)
.
output()
.
/* 0utput of Approximate Nodal Values of u */
void output(double fm[], int nnode, double x[], double y[])
{
int j,m;
FILE *fpout;
.
if((fpout=fopen("outdata.txt","a")) == NULL){
printf("出力データファイルが作成できません．\n");
exit(EXIT_FAILURE);
}
.
.
.
.
.
.
25 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.4 有限要素法プログラム (cnt.)
.
fprintf(fpout, "\nNodal values of u\n");
fprintf(fpout, "
i
u
\n");
for(j=0; j<nnode; j++)
fprintf(fpout, "%4d%12.3e\n", j+1, fm[j]);
fprintf(fpout, "\nMathematica 6 Data\n");
fprintf(fpout, "ListPlot3D[{\n");
for(j=0; j<nnode-1; j++){
fprintf(fpout, "{%f, %f, %f},\n", x[j], y[j], fm[j]);
}
fprintf(fpout, "{%f, %f, %f}\n}, Mesh -> All]\n",
x[j], y[j], fm[j]);
fclose(fpout);
}
.
テキスト [1] からの修正: FILE *fpout; ...
.
.
.
.
.
.
26 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.4 有限要素法プログラム (cnt.)
.
f()
.
/* Function for Free Term */
double f(double x, double y)
{
return(1.0);
}
.
.
.
.
.
.
.
27 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.5 入出力データ
.
indata.txt
.
9 8 5
0 0
0 0.5
0 1
0.5 0
0.5 0.5
0.5 1
1 0
1 0.5
1 1
1 4 5
1 5 2
2 5 6
2 6 3
4 7 8
4 8 5
5 8 9
5 9 6
1
. 2 3 4 7
.
outdata.txt
.
Basic data
nnode =
9
nelmt =
8
nbc
=
5
x,y-coordinates of nodes
i
x(i)
y(i)
1 0.0000 0.0000
2 0.0000 0.5000
3 0.0000 1.0000
4 0.5000 0.0000
5 0.5000 0.5000
6 0.5000 1.0000
7 1.0000 0.0000
8 1.0000 0.5000
. 9 1.0000 1.0000
.
Nodes of elements
i i1 i2 i3
1
1
4
5
2
1
5
2
3
2
5
6
4
2
6
3
5
4
7
8
6
4
8
5
7
5
8
9
5
9
6
. 8
.
.
.
.
.
.
28 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.5 入出力データ (cnt.)
.
Nodes with zero Dirichlet data
1
2
3
4
7
Nodal
i
1
2
3
4
5
6
7
8
. 9
values of u
u
0.000e+000
0.000e+000
0.000e+000
0.000e+000
1.771e-001
2.292e-001
0.000e+000
2.292e-001
3.125e-001
.
Mathematica 6 Data
ListPlot3D[{
{0.000000, 0.000000,
{0.000000, 0.500000,
{0.000000, 1.000000,
{0.500000, 0.000000,
{0.500000, 0.500000,
{0.500000, 1.000000,
{1.000000, 0.000000,
{1.000000, 0.500000,
{1.000000, 1.000000,
},
. Mesh -> All]
.
.
.
0.000000},
0.000000},
0.000000},
0.000000},
0.177083},
0.229167},
0.000000},
0.229167},
0.312500}
.
.
.
29 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.5 入出力データ (cnt.)
Mathematica を立ち上げ，outdata.txt の ListPlot3D[...] を貼り付け，
Shift + Enter により，近似関数を可視化できる．
0.3
1.0
0.2
0.1
0.0
0.5
0.0
0.5
1.0
0.0
.
.
.
.
.
.
30 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.6 演習問題
1.
教材のプログラムを用いて 2 次元 Poisson 問題の数値解を計算し，
Mathematica を用いて結果を図示せよ．ただし，図のような要素分
割を用いよ．基本境界条件と自然境界条件は任意に設定してよい．
.
.
.
.
.
.
31 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.6 演習問題 (cnt.)
2.
授業のウェブページにあるファイル indata 105.txt は，メッシュ作
成プログラムによって，図のようなメッシュを作成し，femfp.c の
入力形式に合わせたデータである．femfp.c で解析し，結果を確認
せよ．
.
.
.
.
.
.
32 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
§7.7 まとめ
2 次元 Poisson 問題を有限要素法で解くためのプログラムの構成につ
いてみてきた．
1.
2.
3.
4.
C 言語で書かれたプログラムは，コンパイル，リンクの手続きによ
り実行可能なプログラムが作成される．
C プログラムの一般形式は関数で構成される．
有限要素法プログラムは，入力，要素係数行列の作成，全体係数行
列の作成，既知項ベクトルの作成，連立 1 次方程式のソルバー，出
力の関数で構成される．
有限要素法プログラムの入力データ形式に合わせて入力ファイルを
作成し，有限要素法プログラムの実行ファイルを走らせて，数値解
が得られることを確認した．
.
.
.
.
.
.
33 / 34
はじめに
有限要素法プログラムの構成
ファイル入出力関数
有限要素法プログラム
入出力データ
演習問題
まとめ
参考文献
参考文献
[1] 菊地文雄.
有限要素法概説 : 理工学における基礎と応用.
サイエンス社, 1980.
.
.
.
.
.
.
34 / 34

Download Report