１.６　有限要素法のプログラム(その１)

これから有限要素法のプログラムを作っていきたいと思います。言語はC++を使いますが、C言語の知識があれば十分理解できると思います。もしC言語の知識がない場合は、先にC言語の簡単な入門書を読んでおくことをお勧めします。
C++はC言語をオブジェクト指向化したもので、クラスがプログラムの構成単位となります。オブジェクト指向については特に説明しませんが、まずこのクラスの例を示します。


class Node
{
    int       ID;
    double    coord[3];          //座標
    double*   value;             //値

public:
    Node()
    {
        value   = NULL;
    }
    virtual ~Node()
    {
        if(value   != NULL) delete[] value;
    }
    void setCoord(int i, double coordi) { coord[i] = coordi; }
    double getCoord(int i) { return coord[i]; }
};

これは節点を定義するためのクラスで最初のキーワード class の後にクラス名 Node とかきます。クラスの定義はそのあとの $\{$ と $\}$ の間にかきます。
この例では三つの変数 ID、coord そして value が定義されています。最初は整数型の変数で「節点番号」、次の倍精度型の配列は座標 $(x,y,z)$ を表し、最後の倍精度のポインター型変数はこの節点が持っている値などを格納するためのものです。これらの変数のことをこのクラスのメンバー変数とよびます。

クラスはこれらの変数だけでなくこれら変数を操作するための関数もメンバー関数として持っています。この例では Node()、^~Node()、setCoord(int i,double coordi)、getCoord(int i)などです。
ここで最初の二つの関数は特別な関数で、Node()はクラス名と同じ名前を持ち関数の型宣言がありません。この関数はこのクラスの変数が作成される時に自動的呼び出される関数で「コンストラクター」といいます。この例ではポインター型変数 value に NULL を代入して未定義になることを避けています。
次の ^~Node はまえに virtual という宣言をされていますがやはり型宣言がありません。これは「デストラクター」とよばれコンストラクターとは逆にこのクラスの型宣言をされた変数が消去されるときに自動的に呼び出されます。この例ではポインター型変数が作成されているときにはそれを消去する処理を行っています。この virtual という宣言は仮想関数であるということですがデストラクターには付けるものだと思っておいてもいいと思います。

さて、これらの関数の定義の前に public というキーワードがありますがこの説明をします。そもそもC言語にはなかったクラスというものを導入したのは、出来るだけ関連した処理はプログラムの中で独立して扱おうということからです。そのためにはこのクラス単位の処理はこの中だけで完結すれば、プログラムは非常に分かりやすくなります。というのもプログラムは大きくなればなるほどいろんな部分が関連して、一つの個所を修正するのにもいろんな場所に影響を与えてしまいます。そこで独立性の高い部品であるクラスと、クラスどうしの関係からプログラムを構築するするということが考えられたのです。

したがってクラスのメンバーである変数や関数はできるだけ外部から直接参照できないことが望ましいのです。そのためにこれらのメンバー変数や関数にどの程度外部からの参照を許すかの基準を設けています。C++ではこの基準を以下の３種類のキーワードで決めています。

private      :　何もしないとこれが適用され、このクラス内の関数からしか参照できません。
public       :　どこからでも参照できます。
protected  :　後から出てくる継承によってこのクラスから派生したクラスからだけ参照できます。

これらのキーワードの後にかかれたメンバーにはこれが適用されます。
この例では座標値を設定したり座標値を外部から参照できる関数が public で定義されています。

このクラスでは関数の実装もこの関数で行っていますが、通常はクラス内部では関数の型宣言だけをしておき、実装は外部で行います。その場合はクラス名を頭につけ、


void Node::setCoord(int i, double coordi)
{
    coord[i] = coordi;
}

などとします。

それではここからプログラムを作っていくことにします。実際に有限要素法のプログラムを作成して、有限要素法の基本的な問題を解けるようにします。
C++のシステムがあればそれを使えば良いのですが、ここでは Windows 環境を使っているものとして無料の MinGW(Minimalist GNU for Windows) がインストールされていることを前提にお話しします。この開発環境は MinGW の公式サイトから入手できます。
まず FEM というホルダーを作成して下さい。そしてその中に次のようなホルダーを作成します。

Include　：クラスを定義したヘッダーファイルを格納するホルダー
Source　 : ソースコードを格納するホルダー

次にメモ帳などのテキストエディターを使ってファイル fem.h をホルダー Include の中に作成し、ファイル fem.cpp を、ホルダー Source の中に作成します。
最初にファイル fem.h の先頭に次のようにヘッダーファイルのインクルード文と定義文をかいてください。



#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <string.h>
#include <sstream>
#include <math.h>

using namespace std;

#define MAX_COLUMN 512
#define ELMAX1  24
#define ELMAX2  24
#define PI 3.14159265358979
#define ZERO_EPS  1.0e-20
#define SMALL_EPS 1.0e-9

次にテキスト読込用クラス Input を定義します。


class Input
{
    //カンマ区切りの文字列の読み込み用クラス
public:
    static int readln(ifstream* fp_in, string strArray[]);
};

このクラスにはメンバー変数はなく、一つのメンバー関数 readln のみを持っています。この関数の実装はファイル fem.cpp で行います。
まず、このファイルに次のインクルード文をかいてからその下にこの関数の実体をかいてください。


#include "../include/fem.h"
#include "../include/solver.h"

int Input::readln(ifstream* fp_in, string strArray[])
{
    //ファイルポインタfp_inから１行読み込みカンマで分割された文字列を返す
    //データの最後はカンマをつける（空文はnum=0）

    int num = 0;

    char str[MAX_COLUMN];
    char buf[MAX_COLUMN];

    while (true)
    {
        fp_in->getline(str, MAX_COLUMN);    //１行読み込み

        int len = strlen(str);              //文字の数
        if(len >= 2)
        {
            //コメント文処理
            if(str[0] == '/' && str[1] == '/')
                continue;
        }
        num = 0;
        int col = 0;
        for (int n=0; n<len+1; n++)
        {
            if(str[n] == ',')
            {
                buf[col] = '\0';
                strArray[num] = buf;
                num++;
                col = 0;
            }
            else
            {
                buf[col] = str[n];
                col++;
                if(str[n] == '\0')
                {
                    strArray[num] = buf;
                    num++;
                    break;
                }
            }
        }
        break;
    }
    return num;
}

次にアドレスクラスをファイル fem.h に定義します。


class Address
{
protected:
    int* address;

public:
    Address(int dim)
    {
        //節点要素（スカラー）  dim = 1
        //節点要素（ベクトル）  dim = 3
        //辺要素                dim = 1
        address = new int[dim];
    }
    virtual ~Address()
    {
        delete[] address;
    }

    void setAddress(int address) { this->address[0] = address; }
    void setAddress(int i, int address) { this->address[i] = address; }
    int getAddress() { return address[0]; }
    int getAddress(int i) { return address[i]; }
};

このクラスは有限要素法の連立方程式を作るときの未知数の番号を定めるもので、スカラー場の場合節点に対応してひとつづつアドレスが与えられますが、拘束条件などがある場合その節点はこの番号から除かれます。またベクトル場の場合節点にベクトルの成分の数だけ番号が付けられます。
このクラスのメンバー関数は全てこの中で実装されていますので、ファイル fem.cpp での記述は必要ありません。

テーブルクラスをファイル fem.h に定義します。


class Table
{
    //テーブルクラス
    int      ID;
    string   title;
    int      type;            // = 1 : 非線形   =  2 : 時間     =  3 : 温度
    int      numtabpt;        //　テーブル点数
    double   *x,*y;
public:
    Table()
    {
        numtabpt = 0;
        x = NULL;
        y = NULL;
    }
    virtual ~Table()
    {
        if(numtabpt > 0)
        {
            delete[] x;
            delete[] y;
        }
    }
    int  getID() { return ID; }
    string getTitle() { return title; }
    int getType() { return type; }
    int getNumtabpt() { return numtabpt; }
    double getX(int i) { return x[i]; }
    double getY(int i) { return y[i]; }
    double getValue(double inx, double& dydx);
    void read(ifstream* fp_in);
    void write(ofstream* fp_out);
};

このクラスは過渡応答解析の時間テーブルや非線形材料の材料特性を表すテーブルを定義するときに使います。メンバー変数には ID やタイプ、テーブル点数と二つの実数を格納する倍精度のポインターを持っています。メンバー関数にはこのテーブルから値を得るための getValue やデータの入出力用のものがあり、実装はファイル fem.cpp で次のように行います。


//テーブルクラス
double Table::getValue(double inx, double& dydx)
{
    //テーブルの範囲外の場合最後の値から外挿した値を使用
    int    n;
    double value = 0.0;

    if(numtabpt <= 0)
    {
        cout << "テーブルエラー   テーブルがセットされていません！" << endl;
        dydx  = 0.0;
        value = 0.0;
        return value;
    }
    if(inx < x[0])
    {
        if(type == 2)    //時間テーブル
        {
            value = 0.0;
            return value;
        }
        cout << "テーブルエラー   " << inx  << " は範囲外の値です！" << endl;
        cout << "                 " << x[0] << " から外挿した値を使用しました！" << endl;
        dydx  = (y[1]-y[0])/(x[1]-x[0]);
        value = y[0] + dydx*(inx-x[0]);
        return value;
    }
    if(x[numtabpt-1] < inx)
    {
        cout << "テーブルエラー   " << inx  << " は範囲外の値です！" << endl;
        cout << "                 " << x[numtabpt-1] << " から外挿した値を使用しました！" << endl;
        dydx  = (y[numtabpt-1]-y[numtabpt-2])/(x[numtabpt-1]-x[numtabpt-2]);
        value = y[numtabpt-1] + dydx*(inx-x[numtabpt-1]);
        return value;
    }
    if(numtabpt == 1)
    {
        if(x[0] > 0.0)
            dydx = y[0]/x[0];
        else
            dydx = 0.0;
        value   = y[0];
    }
    for (n=1; n<numtabpt; n++)
    {
        if((x[n]-x[n-1]) < 1.0e-20)
        {
            cout << "テーブルエラー   Ｘ値 " << x[n] << " は単調増加になっていません！" << endl;
        }
        if(inx <= x[n])
        {
            dydx  = (y[n]-y[n-1])/(x[n]-x[n-1]);
            value = y[n-1] + dydx*(inx-x[n-1]);
            break;
        }
    }
    return value;
}
void Table::read(ifstream* fp_in)
{
    string strArray[MAX_COLUMN];
    char str[MAX_COLUMN];

    //ID、タイプ
    Input::readln(fp_in, strArray);
    strcpy(str,strArray[0].c_str()); ID   = atoi(str);
    strcpy(str,strArray[1].c_str()); type = atoi(str);

    //タイトル
    Input::readln(fp_in, strArray);
    title = strArray[0];

    //テーブル点数
    Input::readln(fp_in, strArray);
    strcpy(str,strArray[0].c_str()); numtabpt = atoi(str);

    if(numtabpt > 0)
    {
        x = new double[numtabpt];
        y = new double[numtabpt];
    }

    for (int i=0; i<numtabpt; i++)
    {
        Input::readln(fp_in, strArray);
        strcpy(str,strArray[0].c_str()); x[i] = atof(str);
        strcpy(str,strArray[1].c_str()); y[i] = atof(str);
    }
}
void Table::write(ofstream* fp_out)
{
    *fp_out << ID << "," << type << "," << endl;
    *fp_out << title << endl;
    *fp_out << numtabpt << "," << endl;
    for (int i=0; i<numtabpt; i++)
        *fp_out << x[i] << "," << y[i] << "," << endl;
}

物性を定義するクラス Material をファイル fem.h にかきます。


class Material
{
    //物性クラス
    int     ID;
    string  title;
    double  value[2];
    Table*  table[2];

public:
    int getID() {return ID; }
    string getTitle() { return title; }
    double getValue(int i) { return value[i]; }
    Table* getTable(int i) { return table[i]; }
    void read(int numtab, Table** tableArray, ifstream* fp_in);
    void write(ofstream* fp_out);
};

線形材料の場合は値をメンバー変数 value に格納します。２個の配列を持っているのは例えば電気伝導率と非透磁率などを指定するためです。非線形材料の場合は Table クラスの変数 table に材料特性テーブルを入れます。
データの入出力用メンバー変数の実装はファイル fem.cppp で行います。


//物性クラス
void Material::read(int numtab, Table** tableArray, ifstream* fp_in)
{
    string strArray[MAX_COLUMN];
    char str[MAX_COLUMN];

    //ID
    Input::readln(fp_in, strArray);
    strcpy(str,strArray[0].c_str()); ID   = atoi(str);

    //タイトル
    Input::readln(fp_in, strArray);
    title = strArray[0];

    //物性値
    Input::readln(fp_in, strArray);
    strcpy(str,strArray[0].c_str()); value[0] = atof(str);
    strcpy(str,strArray[1].c_str()); value[1] = atof(str);

    //非線形テーブル
    int tableID[2];
    Input::readln(fp_in, strArray);
    strcpy(str,strArray[0].c_str()); tableID[0] = atoi(str);
    strcpy(str,strArray[1].c_str()); tableID[1] = atoi(str);
    for (int i=0; i<2; i++)
    {
        table[i] = NULL;
        for (int n=0; n<numtab; n++)
        {
            if(tableArray[n]->getID() == tableID[i])
            {
                table[i] = tableArray[n];
                break;
            }
        }
    }  
}
void Material::write(ofstream* fp_out)
{
    *fp_out << ID << "," << endl;
    *fp_out << title << endl;
    *fp_out << value[0] << "," << value[1] << endl;
}

境界条件のためのクラス Boundary と荷重のためのクラス Load をファイル fem.h に定義します。


class Boundary
{
    //境界条件クラス
    int    ID;
    double value;

public:
    int getID() {return ID; }
    double getValue() { return value; }
    void read(ifstream* fp_in);
    void write(ofstream* fp_out);
};

class Load
{
    //荷重クラス
    int    ID;
    int    type;
    double value[3];
    Table* table;

public:
    int getID() {return ID; }
    int getType() { return type; }
    double getValue(int i) { return value[i]; }
    Table* getTable() { return table; }
    void read(int numtab, Table** tableArray, ifstream* fp_in);
    void write(ofstream* fp_out);
};

境界条件は節点や辺が拘束されている場合に参照されるものでメンバー変数として ID と拘束値 value を持っています。
また荷重クラスは電流密度や磁化などのベクトル値を指定するために３個の配列と時間依存性を表すテーブルを持っています。
これらのクラスの入出力関数はファイル fem.cpp で次のように実装します。


//境界条件クラス
void Boundary::read(ifstream* fp_in)
{
    string strArray[MAX_COLUMN];
    char str[MAX_COLUMN];

    //ID
    Input::readln(fp_in, strArray);
    strcpy(str,strArray[0].c_str()); ID   = atoi(str);

    //境界値
    Input::readln(fp_in, strArray);
    strcpy(str,strArray[0].c_str()); value = atof(str);
}
void Boundary::write(ofstream* fp_out)
{
    *fp_out << ID << "," << endl;
    *fp_out << value << "," << endl;
}

//荷重クラス
void Load::read(int numtab, Table** tableArray, ifstream* fp_in)
{
    string strArray[MAX_COLUMN];
    char str[MAX_COLUMN];

    //ID、タイプ
    Input::readln(fp_in, strArray);
    strcpy(str,strArray[0].c_str()); ID   = atoi(str);
    strcpy(str,strArray[1].c_str()); type = atoi(str);

    //荷重値
    Input::readln(fp_in, strArray);
    strcpy(str,strArray[0].c_str()); value[0] = atof(str);
    strcpy(str,strArray[1].c_str()); value[1] = atof(str);
    strcpy(str,strArray[2].c_str()); value[2] = atof(str);

    //時間テーブル
    int tableID;
    Input::readln(fp_in, strArray);
    strcpy(str,strArray[0].c_str()); tableID = atoi(str);
    table = NULL;
    for (int n=0; n<numtab; n++)
    {
        if(tableArray[n]->getID() == tableID)
        {
            table = tableArray[n];
            break;
        }
    }
}
void Load::write(ofstream* fp_out)
{
    *fp_out << ID << "," << type << "," << endl;
    *fp_out << value[0] << "," << value[1] << "," << value[2] << "," << endl;
}

この節最後に、節点のためのクラス Node をファイル fem.h に定義します。


class ElementList;
class Node
{
    //節点クラス
    int       ID;
    double    coord[3];
    Address*  address;
    Boundary* boundary;
    Load*     load;
    double*   value;

    ElementList* elementlist;

public:
    Node()
    {
        address = NULL;
        boundary = NULL;
        load = NULL;
        value = new double[2];
        for (int i=0; i<2; i++)
            value[i] = 0.0;
        elementlist = NULL;
    }
    virtual ~Node()
    {
        delete[] value;
    }
    int getID() {return ID; }
    double getCoord(int i) { return coord[i]; }
    void setAddress(Address* address) { this->address = address; }
    void setAddress(int i, int address) { this->address->setAddress(i,address); }
    int getAddress()
    {
        if(address != NULL)
            return address->getAddress();
        else
            return -1;
    }
    int getAddress(int i)
    {
        if(address != NULL)
            return address->getAddress(i);
        else
            return -1;
    }
    Boundary* getBoundary() { return boundary; }
    Load* getLoad() { return load; }
    void setValue(int i, double value) { this->value[i] = value; }
    double* getValue() { return value; }
    double getValue(int i) { return value[i]; }
    void setElementlist(ElementList* elementlist) { this->elementlist = elementlist; }
    ElementList* getElementlist() { return elementlist; }
    void read(int numbc, Boundary** boundArray, int numld, Load** loadArray, ifstream* fp_in);
    void write(ofstream* fp_out);
};

最初に class ElementList; とかかれているのはこのクラスのメンバー変数にこのクラスの型が使われているからです。
このクラスはこれより後で定義しているためこの記述がないと未定義の型としてエラーとみなされるからです。
このクラスはメンバー変数として ID と座標値、アドレス、境界条件、荷重を持っています。またこの節点としての値 value をもっています。
メンバー関数である入出力関数の実装はファイル fem.cpp で次のように行います。


//節点クラス
void Node::read(int numbc, Boundary** boundArray, int numld, Load** loadArray, ifstream* fp_in)
{
    string strArray[MAX_COLUMN];
    char str[MAX_COLUMN];

    int boundID,loadID;

    //ID、座標値、荷重条件
    Input::readln(fp_in, strArray);
    strcpy(str,strArray[0].c_str());      ID   = atoi(str);
    strcpy(str,strArray[1].c_str()); boundID   = atoi(str);
    strcpy(str,strArray[2].c_str()); loadID    = atoi(str);
    strcpy(str,strArray[3].c_str()); coord[0]  = atof(str);
    strcpy(str,strArray[4].c_str()); coord[1]  = atof(str);
    strcpy(str,strArray[5].c_str()); coord[2]  = atof(str);
    boundary = NULL;
    for (int n=0; n<numbc; n++)
    {
        if(boundArray[n]->getID() == boundID)
        {
            boundary = boundArray[n];
            break;
        }
    }
    load = NULL;
    for (int n=0; n<numld; n++)
    {
        if(loadArray[n]->getID() == loadID)
        {
            load = loadArray[n];
            break;
        }
    }
}
void Node::write(ofstream* fp_out)
{
    int boundID = 0;
    if(boundary != NULL)
        boundID = boundary->getID();
    int loadID = 0;
    if(load != NULL)
        loadID = load->getID();
    *fp_out << ID << "," << boundID  << "," << loadID << "," << coord[0] << "," << coord[1] << "," << coord[2] << "," << endl;
}

１.６ 有限要素法のプログラム(その１)

１.６　有限要素法のプログラム(その１)