ばねモデルの要素剛性マトリックス、全体剛性マトリックスを作るプログラムを作成してみたの巻

 

1次元のばね要素の剛性マトリックスは、ばねの比例定数をＫとすると

 　K  -K

　-K  K

なので、後はこれを重ね合わせればよいのみ（＾＾）

プログラム例は

public class SpringMain {

    public static void main(String[] args) {

        int nElement = 3; // 要素数
        int nNode = 4; // 節点数
        int nDof = nNode; // 自由度数
        int nSpringElemNode = 2; // ばね要素一つが持つ節点の数
        int[][] elementNode = new int[nElement][nSpringElemNode]; // 要素に属する節点番号

        double[] springK = new double[nElement]; // 各要素のばね定数
        double[][] elemK = new double[2][2]; // 要素剛性マトリックス
        double[][] totalK = new double[nDof][nDof]; // 全体剛性マトリックス

        // 要素が持つ節点の定義
        elementNode[0][0] = 0;
        elementNode[0][1] = 1;

        elementNode[1][0] = 1;
        elementNode[1][1] = 2;

        elementNode[2][0] = 2;
        elementNode[2][1] = 3;

        // ばね定数定義
        springK[0] = 10.0;
        springK[1] = 30.0;
        springK[2] = 20.0;

        for (int i = 0; i < nElement; i++) {

            // 要素剛性マトリックスの作成
            double springKi = springK[i];
            buildElementK(springKi, elemK);

            // 全体剛性マトリックスの作成
            buildTotalK(elementNode, elemK, totalK, i);
        }

        // 全体剛性マトリックスの表示
        for (int i = 0; i < nDof; i++) {
            for (int j = 0; j < nDof; j++) {
                System.out.println("K" + i + j + "=" + totalK[i][j]);
            }
        }

    }

    // 要素剛性マトリックス作成メソッド
    private static void buildElementK(double springK, double[][] elemK) {

        elemK[0][0] = springK;
        elemK[0][1] = -1 * springK;
        elemK[1][0] = -1 * springK;
        elemK[1][1] = springK;
    }

    // 全体剛性マトリックス作成メソッド
    private static void buildTotalK(int[][] elemNode, double[][] elemK,
            double[][] totalK, int N) {

        totalK[elemNode[N][0]][elemNode[N][0]] = totalK[elemNode[N][0]][elemNode[N][0]]
                + elemK[0][0];
        totalK[elemNode[N][0]][elemNode[N][1]] = totalK[elemNode[N][0]][elemNode[N][1]]
                + elemK[0][1];
        totalK[elemNode[N][1]][elemNode[N][0]] = totalK[elemNode[N][1]][elemNode[N][0]]
                + elemK[1][0];
        totalK[elemNode[N][1]][elemNode[N][1]] = totalK[elemNode[N][1]][elemNode[N][1]]
                + elemK[1][1];
    }

}

FEMの簡単な問題設定してみるの巻

 

有限要素法プログラムを作っていくわけだけど、

さすがにいきなり汎用的なのは難しいので、

簡単な問題を解くものから作っていこうと思います。

と言うわけでばね系

 

 

を考えて、これで剛性マトリックスを作ってみましょう。

有限要素法（ＦＥＭ）プログラミングをＪＡＶＡで行ってみることを考える

 

ということで、Ｊａｖａの数値演算ライブラリ

jlapack

を使って有限要素プログラミングが

簡単

にできないか、考えてみた。

ＦＥＭの流れは、基本的に

データ入力

↓

要素剛性マトリックスの作成

↓

全体剛性マトリックスの作成

↓

境界条件自由度処理

↓

連立方程式計算

↓

結果データ処理、出力

なので、こらなければそんなに難しくないはず。

というわけで考えて見ます。

 

 

 

 

 

jblasで行列の足し算と引き算の巻

 

あと、行列の足し算と掛け算ができるようになっておく。

コレにはjlapackではなくて、jblasをつかうとのこと。

jblasは基礎的な演算、jlapackはjblasを組み合わせてできる複雑な演算用とのことらしい。

 

まず足し算から。。。

といっても行列の足し算はないらしい。。。

でも、jlapackでは2次元行列を1次元配列で表現しているので、

ベクトルの足し算ができればよいのである、とのこと。。。

ということで、クラスDapxyを使用

public class DaxpyTest {
    public static void main(String[] args) {
        int N = 4;
        int INCX = 1;
        int INCY = 1;

        double DA = 1.0;
        double[] DX = new double[4];
        double[] DY = new double[4];

        DX[0] = 34.0;
        DX[1] = 25.0;
        DX[2] = 24.0;
        DX[3] = 25.0;
        DY[0] = 25.0;
        DY[1] = 30.0;
        DY[2] = 26.0;
        DY[3] = 23.0;

        Daxpy.daxpy(N, DA, DX, 0, INCX, DY, 0, INCY);

        System.out.println("Answer = ");
        for (int i = 0; i < DY.length; i++)
            System.out.print(DY[i] + " ");
        System.out.println();
    }
}

 

次に掛け算

これもjblasのdgemmを使用。

Ｃマトリックスは使わないが、ダミーとして入ってます。

package test01;

import org.netlib.blas.Dgemm;

public class DgemmTest {
    public static void main(String[] args) {

        int M = 4;
        int N = 4;
        int K = 4;
        double ALPHA = 1.0;

        double[] A = new double[16];
        int LDA = 4;
        double[] B = new double[16];
        int LDB = 4;

        double BETA = 0;
        double[] C = new double[16];
        int LDC = 4;

        A[0] = 2.0;
        A[1] = 1.0;
        A[2] = 5.0;
        A[3] = 2.0;
        A[4] = 1.0;
        A[5] = 4.0;
        A[6] = 3.0;
        A[7] = 2.0;
        A[8] = 5.0;
        A[9] = 3.0;
        A[10] = 1.0;
        A[11] = 3.0;
        A[12] = 2.0;
        A[13] = 2.0;
        A[14] = 3.0;
        A[15] = 2.0;

        B[0] = 34.0;
        B[1] = 25.0;
        B[2] = 24.0;
        B[3] = 25.0;
        B[4] = 25.0;
        B[5] = 30.0;
        B[6] = 26.0;
        B[7] = 23.0;
        B[8] = 24.0;
        B[9] = 26.0;
        B[10] = 44.0;
        B[11] = 25.0;
        B[12] = 25.0;
        B[13] = 23.0;
        B[14] = 25.0;
        B[15] = 21.0;

        C[0] = 34.0;
        C[1] = 25.0;
        C[2] = 24.0;
        C[3] = 25.0;
        C[4] = 25.0;
        C[5] = 30.0;
        C[6] = 26.0;
        C[7] = 23.0;
        C[8] = 24.0;
        C[9] = 26.0;
        C[10] = 44.0;
        C[11] = 25.0;
        C[12] = 25.0;
        C[13] = 23.0;
        C[14] = 25.0;
        C[15] = 21.0;

        Dgemm.dgemm("N", "N", M, N, K, ALPHA, A, 0, LDA, B, 0, LDB, BETA, C, 0,
                LDC);

        System.out.println("Answer = ");
        for (int i = 0; i < C.length; i++)
            System.out.print(C[i] + " ");
        System.out.println();
    }

}

基本的にコレだけ行列計算ができれば、

ＣＡＥ計算プログラムが作れるはず。。。(^^;)

 

 

 

 

 

 

 

 

 

 

jlapackで連立一次方程式を解く

まずはdgesvを使ってみる。
Eclipseの設定でjlapackはインポートしておく

public class DgesvTest {
    public static void main(String[] args) {
        int N = 4;
        int NRHS = 1;
        int LDA = 4;
        int LDB = 4;

        double[] A = new double[16];
        double[] B = new double[4];

        intW info = new intW(0);
        ;

        int[] IPIV = { 0, 0, 0, 0 };

        A[0] = 1.0;
        A[1] = 1.0;
        A[2] = 2.0;
        A[3] = 1.0;
        A[4] = 2.0;
        A[5] = 1.0;
        A[6] = 1.0;
        A[7] = 0.0;
        A[8] = 1.0;
        A[9] = -2.0;
        A[10] = -1.0;
        A[11] = 0.0;
        A[12] = 4.0;
        A[13] = 8.0;
        A[14] = 8.0;
        A[15] = -1.0;
        B[0] = 22.0;
        B[1] = 29.0;
        B[2] = 33.0;
        B[3] = -3.0;

        Dgesv.dgesv(N, NRHS, A, 0, LDA, IPIV, 0, B, 0, LDB, info);

        System.out.println("Answer = ");
        for (int i = 0; i < B.length; i++)
            System.out.print(B[i] + " ");
        System.out.println();
    }
}

行列にイニシャライザを使わないのは0から始まるか1から始まるかわからなくならないようにするため（＾＾）

次にdposvを使ってみる。正定値対称行列用。FEM用にはこっちかな。


public class DposvTest {
public static void main(String[] args) {
int N = 4;
int NRHS = 1;
int LDA = 4;
int LDB = 4;

double[] A = new double[16];
double[] B = new double[4];

intW info = new intW(0);
;

int[] IPIV = { 0, 0, 0, 0 };

A[0] = 34.0;
A[1] = 25.0;
A[2] = 24.0;
A[3] = 25.0;
A[4] = 25.0;
A[5] = 30.0;
A[6] = 26.0;
A[7] = 23.0;
A[8] = 24.0;
A[9] = 26.0;
A[10] = 44.0;
A[11] = 25.0;
A[12] = 25.0;
A[13] = 23.0;
A[14] = 25.0;
A[15] = 21.0;

B[0] = 40.0;
B[1] = 12.0;
B[2] = 22.0;
B[3] = 45.0;

Dposv.dposv("U", N, NRHS, A, 0, LDA, B, 0, LDB, info);
System.out.println("info = " + info + "\n");
System.out.println("Answer = ");
for (int i = 0; i < B.length; i++)
System.out.print(B[i] + " ");
System.out.println();
}
}

あと帯び行列に使えそうなのもあるが、配列の入れ方がよくわからないのでパスしておこう。