Next: 参考文献 Up: LAPACK の dgbsv で連立1次方程式を解く Previous: 5.2 testlapack3.c
(Stokes方程式を有限要素法で解くときのプログラムの書き方は、 公開していなかったっけ… おおむね川原 [2] に従って書いてあります。)
/*
* stokes3.c --- 定常 Stokes 問題の有限要素解を求めるプログラム
* 連立1次方程式を解くために LAPACK を利用
*
* oyabun でのコンパイル
* g77 -I/usr/local/include -O -o stokes3 stokes3.c -llapack -lblas -lmatrix
*/
/* 実行は ./file_name 領域データ 流速の保存ファイル名 圧力の保存ファイル名 */
#include <math.h>
#include <stdio.h>
#include <matrix.h>
#define USELAPACK
void dgbsv_(int *n, int *kl, int *ku, int *nrhs,
double *AB, int *ldab, int *ipiv, double *b, int *ldb, int *info);
int max(int a, int b) { return (a > b) ? a : b; }
int min(int a, int b) { return (a < b) ? a : b; }
/* 一次元配列で計算するのは生で書くと大変なのでマクロで処理 */
#define AB(i,j) ab[(j)*nb+(kl+ku+(i)-(j))]
/* 関数定義*/
/* 全体係数マトリックス全体自由項ベクトルを生成する関数 */
void make_matrix(int T, vector ab, int kl, int ku, int nb,
vector u, double nu, double rho,
int **num, int *cum,
vector fx, vector fy, double hx, double hy)
{
int i,j,k,I,J;
double hx2,hy2,nurho;
matrix A = new_matrix(9,9);
matrix B = new_matrix(9,4);
matrix C = new_matrix(9,4);
matrix D = new_matrix(9,9);
hx2 = hx*hx; hy2 = hy*hy;
nurho = nu * rho;
for (i = 0; i < 9; i++) {
for (j = 0; j < 4; j++)
B[i][j] = C[i][j] = 0;
for (j = 0; j < 9; j++)
A[i][j] = D[i][j] = 0;
}
/* 行列Aの成分の入力*/
for (i=0; i<=3; i++)
A[i][i]= 28* (hx2 + hy2);
A[4][4]=A[6][6]= 112 * hx2 + 64 * hy2;
A[5][5]=A[7][7]= 64 * hx2 + 112 * hy2;
for (i=0;i<=3;i++)
A[i][8]= -16 * ( hx2 + hy2 );
A[4][5]=A[4][7]=A[5][6]=A[6][7]= -16 * ( hx2 + hy2 );
A[4][8]=A[6][8]= -128 * hx2 + 32 * hy2;
A[5][8]=A[7][8]= 32 * hx2 - 128 * hy2;
A[0][5]=A[1][7]=A[2][7]=A[3][5]= 8 * hx2 + 2 * hy2;
A[0][6]=A[1][6]=A[2][4]=A[3][4]= 2 * hx2 + 8 * hy2;
A[0][4]=A[1][4]=A[2][6]=A[3][6]= 14 * hx2 - 32 * hy2;
A[0][7]=A[1][5]=A[2][5]=A[3][7]= -32 * hx2 + 14 * hy2;
A[0][1]=A[2][3]= -7 * hx2 + 4 * hy2;
A[0][3]=A[1][2]= 4 * hx2 - 7 * hy2;
A[0][2]=A[1][3]= -hx2 - hy2;
A[4][6]= 16 * hx2 - 16 * hy2;
A[5][7]= -16 * hx2 + 16 * hy2;
A[8][8]= 256 * hx2 + 256 * hy2;
for (j=0;j<9;j++)
for (i=0;i<9;i++)
A[i][j] = A[j][i];
for (i=0;i<9;i++)
for (j=0;j<9;j++)
A[i][j] = (1.0/(90*hx*hy)) * A[i][j];
/* 行列Bの成分の入力*/
for(i=0;i<9;i++)
for(j=0;j<4;j++)
B[i][j]=0;
B[0][0]=B[3][3]=-25; B[1][1]=B[2][2]=25;
B[1][0]=B[2][3]=5; B[0][1]=B[3][2]=-5;
B[4][0]=B[6][3]=20; B[4][1]=B[6][2]=-20;
B[5][0]=B[5][3]=10; B[5][1]=B[5][2]=50;
B[7][0]=B[7][3]=-50; B[7][1]=B[7][2]=-10;
B[8][0]=B[8][3]=40; B[8][1]=B[8][2]=-40;
for (i=0;i<9;i++)
for (j=0;j<4;j++)
B[i][j] = (hy/180.0) * B[i][j];
/* 行列Cの成分の入力*/
for(i=0;i<9;i++)
for(j=0;j<4;j++)
C[i][j]=0;
C[0][0]=C[1][1]=-25;
C[2][2]=C[3][3]=25;
C[0][3]=C[1][2]=-5;
C[3][0]=C[2][1]=5;
C[4][0]=C[4][1]=-50;
C[4][2]=C[4][3]=-10;
C[5][1]=C[7][0]=20;
C[5][2]=C[7][3]=-20;
C[6][0]=C[6][1]=10;
C[6][2]=C[6][3]=50;
C[8][0]=C[8][1]=40;
C[8][2]=C[8][3]=-40;
for (i=0;i<9;i++)
for (j=0;j<4;j++)
C[i][j] = (hx/180.0) * C[i][j];
/* 行列Dの成分の入力 */
for(i=0;i<=3;i++)
D[i][i]=16;
for(i=4;i<=7;i++)
D[i][i]=64;
for(i=0;i<=3;i++)
D[i][8]=4;
D[4][5]=D[4][7]=D[5][6]=D[6][7]=4;
for(i=4;i<=7;i++)
D[i][8]=32;
D[0][5]=D[0][6]=D[1][6]=D[1][7]=D[2][4]=D[2][7]=D[3][4]=D[3][5]=-2;
D[0][4]=D[0][7]=D[1][4]=D[1][5]=D[2][5]=D[2][6]=D[3][6]=D[3][7]=8;
D[0][1]=D[0][3]=D[1][2]=D[2][3]=-4;
D[0][2]=D[1][3]=1;
D[4][6]=D[5][7]=-16;
D[8][8]=256;
for (j=0;j<9;j++)
for (i=0;i<9;i++)
D[i][j]=D[j][i];
for (i=0;i<9;i++)
for (j=0;j<9;j++)
D[i][j] = (hx*hy/900.0) * D[i][j];
/* 直接剛性法 */
for (k=0;k<T;k++){
for (i=0;i<9;i++){
I = cum[num[k][i]];
for (j=0;j<9;j++) {
J = cum[num[k][j]];
AB(I,J) += nurho * A[i][j];
AB(I+1,J+1) += nurho * A[i][j];
}
for (j=0;j<4;j++){
J = cum[num[k][j]];
AB(I,J+2) -= B[i][j];
AB(I+1,J+2) -= C[i][j];
AB(J+2,I) -= B[i][j];
AB(J+2,I+1) -= C[i][j];
}
for (j=0;j<9;j++) {
u[I] += rho*D[i][j]*fx[num[k][j]];
u[I+1] += rho*D[i][j]*fy[num[k][j]];
}
}
}
free_matrix(A);
free_matrix(B);
free_matrix(C);
free_matrix(D);
}
void change(int n, vector ab, int kl, int ku, int nb,
int k)
{
int i,j;
for (i = max(0,k-ku); i <= min(n-1,k+kl); i++)
AB(i,k) = 0.0;
for (j = max(0,k-kl); j <= min(n-1,k+ku); j++)
AB(k,j) = 0.0;
AB(k,k) = 1.0;
}
void bound(vector ab, int kl, int ku, int nb,
vector u, int *b, vector bx, vector by,
int D, int *cum, int n)
{
int j,i,k;
for (k = 0; k < D; k++) {
j = cum[b[k]];
/* j列を移項 */
for (i = max(0,j-ku); i <= min(n-1,j+kl); i++)
u[i] -= AB(i,j) * bx[k];
u[j] = bx[k];
change(n, ab, kl, ku, nb, j);
/* j+1列を移項 */
j++;
for (i = max(0,j-ku); i <= min(n-1,j+kl); i++)
u[i] -= AB(i,j) * by[k];
u[j] = by[k];
change(n, ab, kl, ku, nb, j);
}
change(n, ab, kl, ku, nb, 2);
u[2] = 0.0;
}
/* main文*/
int main(int argc, char **argv)
{
/**************** 変数宣言***************/
int i,j;
/* 入力ファイル、出力ファイル*/
FILE *f1,*f2,*f3;
/* nelmt:要素数, nnode:総接点数, num_unknown:未知数の総数, nband:半バンド幅 */
int nnode,nelmt,num_unknown,nband,cum_tmp,max_cum,min_cum;
/* nDir: Dirichlet境界上の節点の総数*/
int nDir;
int *cum, **num, *b;
/* nu:粘性定数, rho:密度, hx,hy:x,yについてのメッシュ幅*/
double nu,rho,hx,hy;
/* bx,by:境界値, fx,fy:外力密度 */
vector x,y,bx,by,fx,fy;
/* 全体自由項ベクトル*/
vector u;
/* 全体係数マトリックス*/
int n,ku,kl,nb,ldab,ldb,nrhs,info;
ivector ipvt;
vector ab;
if (argc != 4) {
fprintf(stderr,
"usage: %s <入力データ> <流速データ> <圧力データ>\n", argv[0]);
exit(1);
}
if ((f1 = fopen(argv[1],"r")) == NULL){
fprintf(stderr, "Can't open %s",argv[1]);
exit(1);
}
/* 総要素数, 総節点数を読む */
fscanf(f1, "%d %d %lf %lf", &nelmt, &nnode, &hx, &hy);
printf("要素数,総節点数,hx,hy=\n%d %d %f %f\n", nelmt, nnode, hx, hy);
x = new_vector(nnode);
y = new_vector(nnode);
if (x == NULL || y == NULL) {
fprintf(stderr, "節点座標を記憶するメモリが確保できません。");
exit(1);
}
if ((num = malloc(nelmt * sizeof(void *))) == NULL) {
fprintf(stderr, "要素と節点の対応用のメモリーが足りません\n");
exit(1);
}
for (i=0;i<nelmt;i++) {
num[i] = malloc(sizeof(int) * 9);
if (num[i] == NULL) {
fprintf(stderr, "要素と節点の対応用のメモリーが足りません\n");
exit(1);
}
}
if ((cum = malloc(sizeof(int) * nnode)) == NULL) {
fprintf(stderr, "累積節点番号用のメモリーが足りません \n");
exit(1);
}
/* 節点の座標、要素節点番号対応表、累積節点番号表、総節点数の読み込み */
for (i=0;i<nnode;i++)
fscanf(f1, "%lf %lf", &x[i], &y[i]);
for (i=0;i<nelmt;i++)
for (j=0;j<9;j++)
fscanf(f1, "%d", &num[i][j]);
for (i=0;i<nnode;i++)
fscanf(f1, "%d", &cum[i]); /* i:全体節点番号, cum:累積節点番号 */
fscanf(f1, "%d", &num_unknown);
printf("未知数の個数=%d\n", num_unknown);
if ((u = new_vector(num_unknown)) == NULL) {
fprintf(stderr, "近似解を記憶するメモリーが不足しています。");
exit(1);
}
/* 半バンド幅 nband の評価 */
nband=0;
for (i=0;i<nelmt;i++){
max_cum=0;
min_cum=num_unknown;
for (j=0;j<9;j++) {
cum_tmp=cum[num[i][j]];
if (max_cum < cum_tmp) max_cum = cum_tmp;
if (min_cum > cum_tmp) min_cum = cum_tmp;
}
if (nband < max_cum - min_cum) nband = max_cum - min_cum;
}
printf("半バンド幅=%d\n", nband);
for (i = 0; i < num_unknown; i++)
u[i] = 0;
kl = ku = nband + 2;
nb = 2 * kl + ku + 1;
n = num_unknown;
ab = new_vector(nb * n);
ipvt = new_ivector(n);
if (ab == NULL || ipvt == NULL) {
fprintf(stderr, "係数行列を記憶するメモリーが足りません\n");
exit(1);
}
for (i = 0; i < (nb * n); i++)
ab[i] = 0;
fscanf(f1, "%lf %lf", &nu, &rho);
fscanf(f1, "%d", &nDir);
printf("ν,ρ=\n%f %f\n", nu, rho);
printf("境界上の節点=%d\n", nDir);
if ((b = malloc(sizeof(int) * nDir)) == NULL){
fprintf(stderr, "境界の番号用のメモリーが足りません。\n");
exit(0);
}
bx = new_vector(nDir);
by = new_vector(nDir);
if (bx == NULL || by == NULL) {
fprintf(stderr, "境界値データを読み込むためのメモリーがありません。");
exit(1);
}
for (i = 0; i < nDir; i++)
fscanf(f1, "%d %lf %lf", &b[i], &bx[i], &by[i]); /* 境界値 (u,v) */
fx = new_vector(nnode);
fy = new_vector(nnode);
if (fx == NULL || fy == NULL) {
fprintf(stderr, "外力データを読み込むためのメモリーがありません。");
exit(1);
}
for (i = 0; i < nnode; i++)
fscanf(f1, "%lf %lf", &fx[i], &fy[i]); /* 外力のx,y成分 */
fclose(f1);
/* 全体係数マトリックス、全体自由項ベクトルの生成 */
make_matrix(nelmt,ab,kl,ku,nb,
u,nu,rho,num,cum,fx,fy,hx,hy);
bound(ab,kl,ku,nb,
u,b,bx,by,nDir,cum,num_unknown);
free(num);
free_vector(fx);
free_vector(fy);
free(b); /*free(cum);*/
free_vector(bx); free_vector(by);
/* 連立一次方程式を解く */
nrhs = 1;
ldab = nb;
ldb = n;
dgbsv_(&n, &kl, &ku, &nrhs, ab, &ldab, ipvt, u, &ldb, &info);
if (info == 0)
printf(" successful\n");
else if (info > 0)
printf(" U is singular\n");
else
printf("%d-th argument has illegal value.\n", abs(info));
printf(" return from dgbsv\n");
/* 計算結果をファイルへ出力 */
if ((f2 = fopen(argv[2],"w")) == NULL) {
printf("Can't open %s",argv[2]);
exit(2);
}
if ((f3 = fopen(argv[3],"w")) == NULL) {
printf("Can't open %s",argv[3]);
exit(3);
}
j=0;
for(i=0;i<num_unknown-3;){
fprintf(f2, "%f %f %f %f\n", x[j], y[j], u[i], u[i+1]);
if(cum[j+1]-cum[j]==3){
fprintf(f3, "%f %f %f\n", x[j], y[j], u[i+2]);
i++;
if(x[j+1]!=x[j]){
fprintf(f3,"\n");
}
}
i=i+2;
j++;
}
fprintf(f2, "%f %f %f %f\n", x[nnode-1], y[nnode-1], u[num_unknown-3], u[num_unknown-2]);
fprintf(f3, "%f %f %f\n", x[nnode-1], y[nnode-1], u[num_unknown-1]);
fclose(f2);
fclose(f3);
return 0;
}
桂田 祐史