A. Stokes 方程式を解くプログラム

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

桂田 祐史
2017-10-07