• 検索結果がありません。

第7章 有限要素法のプログラミング

N/A
N/A
Info
ダウンロード
Protected

Academic year: 2021

シェア "第7章 有限要素法のプログラミング"

Copied!
34
0
0

読み込み中.... (全文を見る)

今ダウンロードする ( 34 ページ )

全文

(1)

第

7

章 有限要素法のプログラミング

畔上 秀幸

名古屋大学 情報科学研究科 複雑系科学専攻

(2)

§7.1

はじめに

(目標) 2 次元 Poisson 問題を有限要素法で解くためのプログラムの構成につ いて理解する．

(3)

§7.2

有限要素法プログラムの構成

有限要素法プログラムの具体例をみてみよう．授業のウェブページにあるファ イル femfp.c は，文献 [1] のプログラムを基にして作成された．図のような関数 で構成されている． main( ) input( ) output( ) assem( ) ecm( ) f( ) solve( ) gs{solve( )

(4)

§7.3

ファイル入出力関数

fopen()

#include <stdio.h>

FILE *fopen(char *fname, char *mode);

• fname のファイルを mode (下表) にしたがって開く． • 返却値は，成功したときファイルの番号，失敗したとき NULL となる． mode 作用 mode 作用 r テキストの読み込み rb バイナリの読み込み r+ テキストの読み書き rb+ バイナリファイルの読み書き w テキストの書き込み wb バイナリの書き込み w+ テキストの読み書き wb+ バイナリの読み書き a テキストの追加書き込み ab バイナリの追加書き込み a テキストの読み書き ab+ バイナリの読み書き

(5)

fclose()

#include <stdio.h> int fclose(FILE *fp); • fopen で開いたファイル fp を閉じる． • 返却値は，成功したとき 0，失敗したとき EOF となる．

fprintf()

include <stdio.h>

int fprintf(FILE *fp, char *fmt,...);

• ファイル fp に対して，書式付き文字列 fmt を書き込む．

(6)

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 となる．

(7)

§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

(8)

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);

(9)

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

(10)

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); }

(11)

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

(12)

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"); }

(13)

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, ...

(14)

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; }

(15)

/* 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]; } } }

(16)

/* 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; } } }

(17)

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;

(18)

/* 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; }

(19)

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]; } }

(20)

/* 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];

} }

(21)

文献 [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; }

(22)

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]; }

(23)

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]); }

(24)

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]; }

(25)

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); }

(26)

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); }

(27)

f()

/* Function for Free Term */ double f(double x, double y) {

return(1.0); }

(28)

§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 8 5 9 6

(29)

Nodes with zero Dirichlet data 1 2 3 4 7 Nodal values of u i u 1 0.000e+000 2 0.000e+000 3 0.000e+000 4 0.000e+000 5 1.771e-001 6 2.292e-001 7 0.000e+000 8 2.292e-001 9 3.125e-001 Mathematica 6 Data ListPlot3D[{ {0.000000, 0.000000, 0.000000}, {0.000000, 0.500000, 0.000000}, {0.000000, 1.000000, 0.000000}, {0.500000, 0.000000, 0.000000}, {0.500000, 0.500000, 0.177083}, {0.500000, 1.000000, 0.229167}, {1.000000, 0.000000, 0.000000}, {1.000000, 0.500000, 0.229167}, {1.000000, 1.000000, 0.312500} }, Mesh -> All]

(30)

Mathematica を立ち上げ，outdata.txt の ListPlot3D[...] を貼り付け，Shift + Enter により，近似関数を可視化できる．

(31)

§7.6

演習問題

1 教材のプログラムを用いて 2 次元 Poisson 問題の数値解を計算し， Mathematica を用いて結果を図示せよ．ただし，図のような要素分割を用 いよ．基本境界条件と自然境界条件は任意に設定してよい．

(32)

2 授業のウェブページにあるファイル indata 105.txt は，メッシュ作成プログ ラムによって，図のようなメッシュを作成し，femfp.c の入力形式に合わせ たデータである．femfp.c で解析し，結果を確認せよ．

(33)

§7.7

まとめ

2 次元 Poisson 問題を有限要素法で解くためのプログラムの構成についてみて きた． 1 C 言語で書かれたプログラムは，コンパイル，リンクの手続きにより実行可 能なプログラムが作成される． 2 C プログラムの一般形式は関数で構成される． 3 有限要素法プログラムは，入力，要素係数行列の作成，全体係数行列の作 成，既知項ベクトルの作成，連立 1 次方程式のソルバー，出力の関数で構成 される． 4 有限要素法プログラムの入力データ形式に合わせて入力ファイルを作成し， 有限要素法プログラムの実行ファイルを走らせて，数値解が得られることを 確認した．

(34)

参考文献

[1] 菊地文雄.

有限要素法概説: 理工学における基礎と応用.

参照

今ダウンロードする ( PDF - 34 ページ - 526.66 KB )

関連したドキュメント

CSV ファイルをメモ帳で確認 CSV ファイルを確認 編集するときは テキストエディタで確認するとよいと聞きました CSV ファイルを確認 編集するには テキストエディタのご使用をお奨めします テキストエディタを使用すると 不要な タブ 文字が含まれていないか 余計な カンマ がないか JAN コ

【CSV ファイルをメモ帳で確認】 CSV ファイルを確認・編集するときは、テキストエディタで確認するとよいと聞きました。

マルチワーカー取扱説明書 Rev. 1.00A 2016/6/24 本製品をご利用になる前に ここに記載された注意事項をよく読み必ずお守りください 表示マークについて人が死亡または重傷を負う可能性が高いと想定される内容です 人が死亡または重傷を負う可能性が想定される内容です 人が障害を負う可能性また

Clock Mode Error 動作周波数エラーが発生しました。.

Microsoft PowerPoint - 0ßLª°’¬ÇŽHP² (.pptx

ダウンロードしたファイルを 解凍して自動作成ツール (StartPro2018.exe) を起動します。.

充 実 の 課 外 活 動 で 学 生 生 活 ボ ラ ン テ ィ ア

キャンパスの軸線とな るよう設計した。時計台 は永きにわたり図書館 として使 用され、学 生 の勉学の場となってい たが、9 7 年の新 大

砂川 元 沖縄県

学校の PC などにソフトのインストールを禁じていることがある　そのため絵本を内蔵した iPad

罪 法条競合と包括 論

となってしまうが故に︑

令和４年度版　みんなで読める税の本 令和４年度版　みんなで読める税の本

 ・ ナンバープレートを破損、紛失したとき   ・ 住所、氏名、定置場等に変更があったとき  ・

^#：ArduinoはArduino IDEやGUIとの`$にabすること があります

QRされた .ino ファイルを Arduino に‚き1む ことで、 GUI |}した ƒ+どおりに Arduino を/‡((スタンドアローン})させるこ とができます。. 1)