processingで1次元有限要素解析

所用で有限要素法について勉強することになり、無謀にも使い慣れているProcessingでプログラムを組もうと思いました。その結果できたのが、下のプログラムになります。左端にほぼ任意の値を、右端にほぼ任意の傾きを設定して解くことができました。

連立1次方程式を解く際には、ApacheのCommons Mathを使いました。下のプログラムの場合、fem1dフォルダにcodeというフォルダを作り、さらにその中にcommons-math3-3.x.jarファイルを入れて実行する必要があります。

fem1d.pde

//1次元ポアソン問題[-du(x)/dx=q (MinX<x<MaxX),u(MinX)=0,du(MaxX)/dx=0]を、有限要素法で解く
import org.apache.commons.math3.linear.*; //行列計算ライブラリ


void Boundary(int NodeNum, float Dirichlet, float Neumann) {
  //固定境界条件
  if (Dirichlet!=9999) { //Dirichlet=9999の時は処理しない
    for (int n=0; n<N_nod; n++) {
      B_glo[n][0] -= Dirichlet*A_glo[n][NodeNum]; //定数ベクトルの全成分に移項
      A_glo[NodeNum][n] = 0.0; //行を全て0にする
      A_glo[n][NodeNum] = 0.0; //列を全て0にする
    }
    A_glo[NodeNum][NodeNum] = 1.0; //固定する節点番号のみ1にする
    B_glo[NodeNum][0] = Dirichlet; //任意の値で固定
  }

  //自然境界条件
  if (Neumann!=9999) { //Neumann=9999の時は処理しない
    B_glo[NodeNum][0] += Neumann; //任意の傾きで固定
  }
}


float s; //描画の倍率
double Dirichlet=1;

int N_nod=11, N_ele=N_nod-1; //節点数、要素数
float G_nodX[] = new float[N_nod]; //Global節点のx座標
int Ele_nod[][] = new int[N_ele][2]; //各要素のGlobal節点番号
float L_nodX[][] = new float[N_ele][2]; //各要素のLocal節点のx座標

float diff_d=1.0, cons_f=1.0; //係数D、定数項b
float A_ele[][][] = new float[N_ele][2][2]; //各要素の要素行列
float B_ele[][] = new float[N_ele][2]; //各要素の定数ベクトル
double A_glo[][] = new double[N_nod][N_nod]; //全体行列
double B_glo[][] = new double[N_nod][1]; //全体の定数ベクトル

RealMatrix matA; //全体行列A、定数ベクトルB
RealMatrix vecB; //全体行列A、定数ベクトルB
EigenDecomposition detA; //Aの行列式
RealMatrix matAi, vecU; //Aの逆行列、未知数ベクトルU
float U_glo[] = new float[N_nod]; //全体の未知数ベクトル


void setup() {
  size(600, 500);
  translate(0, height/2);
  textAlign(CENTER);
  background(127);
  s = height/150;


  ////////// データの前処理 //////////
  //離散化したx座標を定義(0～10)
  for (int n=0; n<N_nod; n++) {
    G_nodX[n] = 10*n/(N_nod-1);  //計算領域を等分割
  }

  //各要素のGlobal節点番号
  for (int e=0; e<N_ele; e++) {
    Ele_nod[e][0] = e;
    Ele_nod[e][1] = e+1;
  }

  //各要素のLocal節点座標
  for (int e=0; e<N_ele; e++) {
    for (int i=0; i<2; i++) {
      L_nodX[e][i] = G_nodX[ Ele_nod[e][i] ];
    }
  }


  ////////// 要素行列を計算し、全体行列を作る //////////
  for (int i=0; i<N_nod; i++) { //全体行列を0で初期化
    for (int j=0; j<N_nod; j++) {
      A_glo[i][j] = 0.0;
    }
    B_glo[i][0] = 0.0;
  }
  
  for (int e=0; e<N_ele; e++) { //要素行列の各成分を全体行列に足しこむ
    for (int i=0; i<2; i++) {
      for (int j=0; j<2; j++) { //各要素のGlobal節点と、全体行列の行と列が対応
        A_glo[ Ele_nod[e][i] ][ Ele_nod[e][j] ] += pow(-1, i+1)*pow(-1, j+1) / (L_nodX[e][1] -L_nodX[e][0]);
      }
      B_glo[ Ele_nod[e][i] ][0] += cons_f*(L_nodX[e][1]-L_nodX[e][0])/2;
    }
  }

  //全体行列の確認
  println("Global matrix (constructed)"); 
  for (int i=0; i<N_nod; i++) {
    for (int j=0; j<N_nod; j++) {
      print(nfs((float)A_glo[i][j], 2, 2) +" ");
    }
    println(";  " +nfs((float)B_glo[i][0], 2, 2));
  }
  

  ////////// 境界条件処理 //////////
  Boundary(0, -10, 9999); //左端はディリクレ境界
  Boundary(N_nod-1, 9999, -5); //右端はノイマン境界

  //全体行列の確認
  println("Global matrix"); 
  for (int i=0; i<N_nod; i++) {
    for (int j=0; j<N_nod; j++) {
      print(nfs((float)A_glo[i][j], 2, 2) +" ");
    }
    println("; " +nfs((float)B_glo[i][0], 2, 2));
  }


  ////////// 連立方程式Au=bを解く //////////
  //2次元配列を行列として設定
  matA = new Array2DRowRealMatrix(A_glo);
  vecB = new Array2DRowRealMatrix(B_glo);

  //行列式
  detA = new EigenDecomposition(matA);
  println("detA = " +detA.getDeterminant());

  //行列式が非ゼロ(正則)の時のみ、以降の計算が可能
  print(detA.getDeterminant());

  matAi = MatrixUtils.blockInverse(matA, 0); //逆行列
  vecU = matAi.multiply(vecB); //Au=b → u=Ai*b

  //求めた未知数ベクトルの確認
  println("Unkown vector U = ");
  for (int i=0; i<N_nod; i++) {
    U_glo[i] = (float)vecU.getEntry(i, 0); //求めた未知数を新しい配列に代入
    print(nfs(U_glo[i], 3, 2) +", ");
  }


  ////////// 計算結果の可視化 //////////
  //基準線
  stroke(255, 0, 0);
  strokeWeight(1);
  line(0, 0, width, 0); 

  //節点間の線分
  stroke(255);
  strokeWeight(1);
  for (int i=0; i<N_nod-1; i++) {
    float u1 = U_glo[i], u2 = U_glo[i+1]; //配列Uの値を取得
    line(i*width/(N_nod-1), -u1*s, (i+1)*width/(N_nod-1), -u2*s);
  }

  //節点と番号
  stroke(0);
  strokeWeight(3);
  fill(255);
  for (int i=0; i<N_nod; i++) {
    point(i*width/(N_nod-1), -U_glo[i]*s);
    text(i, i*width/(N_nod-1), -U_glo[i]*s -height/60);
  }
}

えー、見てお分かりかもしれませんが、かなり無理をしてしまいました。
こんなだったらPython使った方が簡単じゃんと思いました。
有限要素法の解説とかも書こうかと思いましたが、プログラムがお粗末なので、一旦出直したいと思います……