ここでファイル fem.h に有限要素法用のクラス Fem を次のように定義します。
class Fem
{
protected:
int nstep,istep;
Table* timetable;
double time,dt;
double theta;
bool lnflg; //線形=true、非線形=false
int fieldType; // = 1 : 勾配場
// = 2 : 回転場
// = 3 : テンソル場
int nfree; //自由度
int* num_column; //係数行列の行あたりの非ゼロ数
int** index_matrix; //インデックス行列
double** matrix; //係数行列
double* vector; //右辺列ベクトル
Analysis* analysis;
ofstream* fp_chk; //チェックライト
public:
Fem(Analysis* analysis, bool lnflg, int fieldType, double eps, ofstream* fp_chk)
{
istep = 1;
timetable = NULL;
time = 0.0;
dt = 0.0;
theta = 2.0/3.0;
num_column = NULL;
index_matrix = NULL;
matrix = NULL;
vector = NULL;
this->fp_chk = fp_chk;
}
virtual ~Fem()
{
delete[] num_column;
for (int n=0; n<nfree; n++)
delete[] index_matrix[n];
delete[] index_matrix;
}
void setNstep(int nstep) { this->nstep = nstep; }
void setTheta(double theta) { this->theta = theta; }
int getNstep();
void setTime(int istep);
virtual void exec() = 0;
virtual void setAddress() {};
virtual void setIndexmatrix() {};
virtual void assem() {};
virtual void output(ofstream* fp_out) = 0;
virtual void renew() {};
};
ここでこのクラスのメンバー変数はこのクラスを継承したクラスから参照できるように全てアクセス制限は protected としています。
変数 lnflg は線形解析か非線形解析かのフラッグ、fieldType は扱う場がポテンシャルの勾配場か回転場かを示すフラッグなどです。
nfree 以下 vector までは連立方程式のための変数で、analysis は解析用クラスの変数です。
このクラスのメンバー関数で、exec は有限要素法の実行用関数、setAddress、setIndexmatrix は連立方程式用配列の作成用です。
また assem は要素行列から全体の係数行列を組み立てるアッセンブル用関数、output はアウトプット用、renew は次の時間ステップのための準備を行う関数です。これらの関数は全て仮想関数で、このクラスから派生したクラスで具体的な形が与えられます。この中で exec と output は純粋仮想関数としており必ず派生クラスで実装する必要があります。
有限要素法は変数を節点にもたせるか辺にもたせるかで節点要素法と辺要素法に分かれることを説明しました。また節点要素法の場合も、節点にポテンシャルなどのスカラー値を持たせる場合と、ベクトルポテンシャルや変位などのベクトル値を持たせる場合があります。
これらの有限要素法を扱うには同じ有限要素法といってもかなり扱いが異なりますので、それぞでの有限要素法クラスをもたせた方が便利です。
そこで、スカラーを節点にもつクラス NodalScalarfem、ベクトルを節点にもつクラス NodalVectorfem、そして辺にベクトルの射影値を持たせるクラス Edgefem を次のように定義します。
class NodalScalarfem : public Fem
{
public:
NodalScalarfem(Analysis* analysis, bool lnflg, int fieldType, double eps,
ofstream* fp_chk) : Fem(analysis,lnflg,fieldType,eps,fp_chk)
{
this->analysis = analysis;
this->lnflg = lnflg;
this->fieldType = fieldType;
this->nstep = nstep;
this->fp_chk = fp_chk;
}
~NodalScalarfem()
{
}
void exec();
void setAddress();
void setIndexmatrix();
void assem();
void output(ofstream* fp_out);
void renew();
};
class NodalVectorfem : public Fem
{
protected:
int dim; //ベクトル場の次元
public:
NodalVectorfem(Analysis* analysis, bool lnflg, int fieldType, double eps,
ofstream* fp_chk) : Fem(analysis,lnflg,fieldType,eps,fp_chk)
{
this->analysis = analysis;
this->lnflg = lnflg;
this->fieldType = fieldType;
this->nstep = nstep;
this->fp_chk = fp_chk;
}
~NodalVectorfem()
{
}
void exec();
void setAddress();
void setIndexmatrix();
void assem();
void output(ofstream* fp_out);
void renew();
};
class Edgefem : public Fem
{
public:
Edgefem(Analysis* analysis, bool lnflg, int fieldType, double eps,
ofstream* fp_chk) : Fem(analysis,lnflg,fieldType,eps,fp_chk)
{
this->analysis = analysis;
this->lnflg = lnflg;
this->fieldType = fieldType;
this->nstep = nstep;
this->fp_chk = fp_chk;
}
~Edgefem()
{
}
void exec();
void setAddress();
void setIndexmatrix();
void assem();
void output(ofstream* fp_out);
void renew();
};
これらのメンバー関数の実装をファイル fem.cpp にかきます。ここではクラス NodalScalarfem のみの実装を行います。
void NodalScalarfem::exec()
{
*fp_chk << endl;
*fp_chk << "時間ステップ " << istep << " 時刻 " << time << " ステップ幅 " << dt << endl;
*fp_chk << endl;
bool cvflg = lnflg; //収束フラッグ
for (int newton=1; newton<=analysis->getNewtonmax(); newton++)
{
//係数行列の作成
assem();
//拘束条件の適用
//連立方程式の解法
int iccg_max = analysis->getICCGmax();
double iccg_eps = analysis->getICCGeps();
Iccg* iccg = new Iccg(nfree,num_column,index_matrix,matrix,vector,iccg_max,iccg_eps,fp_chk);
iccg->solver();
double* x = new double[nfree];
for (int n=0; n<nfree; n++)
x[n] = iccg->getX(n);
//連立方程式用の配列のクリア
delete iccg;
//求まった未知数の節点値へのセット
double v2 = 0.0;
double dv2 = 0.0;
for (int n=0; n<analysis->getNumnp(); n++)
{
Node* node = analysis->getNode(n);
double* value = node->getValue();
int address = node->getAddress();
if(address >= 0)
value[0] += x[address];
else
value[0] = 0.0;
v2 += value[0]*value[0];
if(address >= 0)
dv2 += x[address]*x[address];
}
delete[] x;
if(lnflg) break;
//収束判定
double newton_eps = analysis->getNewtoneps();
double v = sqrt(v2);
double dv = sqrt(dv2);
double rat;
if(v > ZERO_EPS)
rat = dv/v;
else if(dv > ZERO_EPS)
rat = 1.0;
else
rat = 0.0;
*fp_chk << " IT. NO = " << newton << " DX/X = " << rat << endl;
cout << " IT. NO = " << newton << " DX/X = " << rat << endl;
if(rat < newton_eps)
{
//収束
cvflg = true;
break;
}
else
{
cvflg = false;
}
}
}
int Fem::getNstep()
{
//時間テーブルの設定
timetable = analysis->getTimetable();
int nstep;
if(timetable == NULL)
nstep = 1; //静解析
else
nstep = timetable->getNumtabpt(); //動解析
return nstep;
}
void Fem::setTime(int istep)
{
//time、dtの設定
if(timetable == NULL)
{
time = 0.0;
dt = 0.0;
}
else
{
time = timetable->getX(istep-1);
dt = time;
dt = timetable->getX(istep-1) - timetable->getX(istep-2);
}
}
void NodalScalarfem::setAddress()
{
int numnp = analysis->getNumnp();
//節点のアドレスと自由度数nfreeを求める
nfree = 0;
for (int n=0; n<numnp; n++)
{
Node* node = analysis->getNode(n);
if(node->getBoundary() == NULL)
{
Address* address = new Address(1);
address->setAddress(nfree);
node->setAddress(address);
nfree++;
}
}
}
void NodalScalarfem::setIndexmatrix()
{
int numnp = analysis->getNumnp();
int numel = analysis->getNumel();
//節点を持つ要素リストの作成
for (int n=0; n<numel; n++)
{
Element* element = analysis->getElement(n);
for (int i=0; i<element->getNumnode(); i++)
{
Node* node = element->getNode(i);
if(node->getElementlist() == NULL)
{
ElementList* elementlist = new ElementList(element);
node->setElementlist(elementlist);
}
else
{
node->getElementlist()->add(element);
}
}
}
//num_columnの配列確保
num_column = new int[nfree];
index_matrix = new int*[nfree];
for (int n=0; n<numnp; n++)
{
NodeList* nodelist = NULL;
Node* nodei = analysis->getNode(n);
int ii = nodei->getAddress();
if(ii < 0) continue;
ElementList* list = nodei->getElementlist();
while (list != NULL)
{
Element* element = list->getElement();
for (int j=0; j<element->getNumnode(); j++)
{
Node* nodej = element->getNode(j);
int jj = nodej->getAddress();
if(jj < 0) continue;
if(jj <= ii)
{
if(nodelist == NULL)
nodelist = new NodeList(nodej);
else
nodelist->add(nodej);
}
}
list = list->getNext();
}
NodeList* p = nodelist;
num_column[ii] = 0;
while (p != NULL)
{
num_column[ii]++;
p = p->getNext();
}
//index_matrixの配列確保
index_matrix[ii] = new int[num_column[ii]];
p = nodelist;
int jj = 0;
while (p != NULL)
{
index_matrix[ii][jj] = p->getNode()->getAddress() + 1;
jj++;
p = p->getNext();
}
//ソート
if(num_column[ii] == 0) continue;
int* index = new int[num_column[ii]];
for (jj=0; jj<num_column[ii]; jj++)
index[jj] = jj;
IntSort* st = new IntSort(num_column[ii]); // ソート用オブジェクトの作成
st->setData(index,index_matrix[ii]);
st->quick_sort(0,num_column[ii]-1);
st->getIndex(index);
st->getData(index_matrix[ii]);
delete st;
delete[] index;
delete nodelist;
}
//節点の要素リストのクリア
for (int n=0; n<numnp; n++)
{
Node* node = analysis->getNode(n);
delete node->getElementlist();
node->setElementlist(NULL);
}
if(check)
{
//check write
*fp_chk << endl;
*fp_chk << "インデックス行列" << endl;
for (int n=0; n<nfree; n++)
{
*fp_chk << n << " " << num_column[n] << " ";
for (int i=0; i<num_column[n]; i++)
*fp_chk << " " << index_matrix[n][i];
*fp_chk << endl;
}
}
}
void NodalScalarfem::assem()
{
//係数行列のための配列確保
matrix = new double*[nfree];
vector = new double[nfree];
for (int n=0; n<nfree; n++)
{
matrix[n] = new double[num_column[n]];
for (int i=0; i<num_column[n]; i++)
matrix[n][i] = 0.0;
vector[n] = 0.0;
}
int numel = analysis->getNumel();
int anal_type = analysis->getAnaltype();
for (int n=0; n<numel; n++)
{
if(n%1000 == 0)
cout << " N = " << (n+1) << " / " << numel << endl;
Element* element = analysis->getElement(n);
int numnode = element->getNumnode();
//要素行列
double elmatrix[3][ELMAX1][ELMAX2];
// double** C = elmatrix[0];
// double** K = elmatrix[1];
// double** Ks = elmatrix[2];
double F[ELMAX1];
element->elMatrix(elmatrix,F,anal_type,lnflg,time,dt,theta);
double K1[ELMAX1][ELMAX2];
double K2[ELMAX1][ELMAX2];
double K3[ELMAX1][ELMAX2];
if(dt > ZERO_EPS)
{
for (int i=0; i<numnode; i++)
{
for (int j=0; j<numnode; j++)
{
K1[i][j] = elmatrix[0][i][j]/dt - (1.0-theta)*elmatrix[1][i][j];
K2[i][j] = elmatrix[0][i][j]/dt + theta*elmatrix[1][i][j];
K3[i][j] = elmatrix[0][i][j]/dt + theta*elmatrix[2][i][j];
}
}
}
else
{
for (int i=0; i<numnode; i++)
{
for (int j=0; j<numnode; j++)
{
K1[i][j] = -(1.0-theta)*elmatrix[1][i][j];
K2[i][j] = theta*elmatrix[1][i][j];
K3[i][j] = theta*elmatrix[2][i][j];
}
}
}
for (int i=0; i<numnode; i++)
{
int ii = element->getNode(i)->getAddress();
if(ii < 0) continue;
for (int j=0; j<numnode; j++)
{
int jj = element->getNode(j)->getAddress();
if(jj < 0) continue;
if(jj <= ii)
{
for (int k=0; k<num_column[ii]; k++)
{
if(jj+1 == index_matrix[ii][k])
{
if(lnflg)
matrix[ii][k] += K2[i][j];
else
matrix[ii][k] += K3[i][j];
break;
}
}
}
}
vector[ii] += F[i];
if(lnflg)
{
for (int j=0; j<numnode; j++)
{
int jj = element->getNode(j)->getAddress();
if(jj < 0) continue;
vector[ii] += K1[i][j]*element->getNode(j)->getValue(1);
}
}
else
{
for (int j=0; j<numnode; j++)
{
int jj = element->getNode(j)->getAddress();
if(jj < 0) continue;
vector[ii] += K1[i][j]*element->getNode(j)->getValue(1)
- K2[i][j]*element->getNode(j)->getValue(0);
}
}
}
}
if(check)
{
///check write
for (int n=0; n<nfree; n++)
{
*fp_chk << "matrix[" << n << "] = " << num_column[n] << " ";
for (int i=0; i<num_column[n]; i++)
*fp_chk << matrix[n][i] << " ";
*fp_chk << endl;
}
for (int n=0; n<nfree; n++)
{
*fp_chk << "vector[" << n << "] = " << vector[n] << endl;
}
}
}
void NodalScalarfem::output(ofstream* fp_out)
{
*fp_out << endl;
*fp_out << "時間ステップ " << istep << " 時刻 " << time << " ステップ幅 " << dt << endl;
*fp_out << endl;
*fp_out << "-------------- ポテンシャル ------------------------" << endl;
*fp_out << endl;
for (int n=0; n<analysis->getNumnp(); n++)
{
Node* node = analysis->getNode(n);
*fp_out << " " << node->getID() << " " << node->getValue(0) << endl;
}
*fp_out << endl;
*fp_out << "---------------- 勾配ベクトル -----------------------" << endl;
*fp_out << endl;
for (int n=0; n<analysis->getNumel(); n++)
{
Element* element = analysis->getElement(n);
double h[8],d[8][3];
double rst[3];
element->getRST(rst);
element->deriv(rst[0], rst[1], rst[2], h, d);
double gradient[2][3];
element->grad(h, d, gradient);
*fp_out << " " << element->getID() << " (" << -gradient[0][0] << ", " << -gradient[0][1] << ")" << endl;
}
}
void NodalScalarfem::renew()
{
for (int n=0; n<analysis->getNumnp(); n++)
{
Node* node = analysis->getNode(n);
node->setValue(1,node->getValue(0)); //未知数を旧値にセット
node->setValue(0,0.0); //未知数をクリア
}
}