この節では要素を扱うクラス Element について説明します。
有限要素法ではいろんな種類の要素を扱います。例えば２次元では三角形要素や四角形要素、３次元では四面体要素や六面体要素などがあります。
これらの要素を統一的に扱うためにまず全てに共通な Element クラスを定義し、この派生クラスとしてそれぞれの要素のクラスを定義します。
ここでは派生クラスとして、２次元の三角形要素 element23、四角形要素 element24、３次元の四面体要素 element34、六面体要素 element38 を定義します。
メンバー変数としては ID、タイプ、この要素を構成する節点を表す Node 型のポインター、物性を表す Material 型のポインター、荷重を表す Load 型ポインターです。メンバー関数としては、ベクトル演算でよく使われる内積用関数 inner と外積用関数 outer をもちこれらは static の宣言をしているのでこのクラスの変数を作らなくても Element::inner(v1,v2) などとして外部から使うことが出来ます。
また、全ての要素で共通に使う変数 nint、xg、wgt は数値積分のための定数で static として宣言されています。
入出力用関数と、形状関数、微分演算をとる関数、要素行列を作成するためのメンバー関数は仮想関数として頭に virtual をつけて定義します。
以下のコードをファイル fem.h に入力してください。
class Element
{
//要素クラス
protected:
int ID;
int type; // = 23:三角形要素 = 24:四角形要素 = 34:四面体要素 = 38:六面体要素
Node** node;
Material* material;
Load* load;
public:
Element(int ID)
{
this->ID = ID;
material = NULL;
load = NULL;
}
//ベクトル演算
static double inner(double* v1, double* v2)
{
//内積
return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
}
static void outer(double* v1, double* v2, double* v3)
{
//外積
v3[0] = v1[1]*v2[2] - v1[2]*v2[1];
v3[1] = v1[2]*v2[0] - v1[0]*v2[2];
v3[2] = v1[0]*v2[1] - v1[1]*v2[0];
}
//数値積分用定数
static int nint;
static double xg[4][4];
static double wgt[4][4];
static Element* getElement(int ID, int type);
int getID() { return ID; }
int getType() { return type; }
virtual Node* getNode(int i) = 0;
virtual void setMaterial(Material* material) = 0;
virtual Material* getMaterial() = 0;
virtual Load* getLoad() = 0;
virtual int getNumnode() = 0;
virtual void read(string strArray[], int numat, Material** matArray, int numld, Load** loadArray, int numnp, Node** nodeArray, ifstream* fp_in) = 0;
virtual void read(int numnode, string strArray[], int numat, Material** matArray, int numld, Load** loadArray, int numnp, Node** nodeArray, ifstream* fp_in);
void write(int numnode, ofstream* fp_out);
//形状関数
virtual void getRST(double rst[]) = 0;
virtual double deriv(double r, double s, double t, double h[], double d[][3]) = 0;
//微分演算
virtual void grad(double h[], double d[][3], double value[][3]) = 0;
//要素行列
virtual void elMatrix(double elmatrix[][ELMAX1][ELMAX2], double elvector[], int anal_type, int lnflg, double time, double dt, double theta) = 0;
};
class Element23 : public Element
{
//３節点要素（３角形要素）
public:
Element23(int ID) : Element(ID)
{
node = new Node*[3];
for (int i=0; i<3; i++)
node[i] = NULL;
}
~Element23()
{
delete[] node;
}
Node* getNode(int i) { return node[i]; }
void setMaterial(Material* material) { this->material = material; }
Material* getMaterial() { return material; }
Load* getLoad() { return load; }
int getNumnode() { return 3; }
void read(string strArray[], int numat, Material** matArray, int numld, Load** loadArray, int numnp, Node** nodeArray, ifstream* fp_in);
void write(ofstream* fp_out);
//形状関数
void getRST(double rst[]);
double deriv(double r, double s, double t, double h[], double d[][3]);
//微分演算
void grad(double h[], double d[][3], double value[][3]);
//要素行列
void elMatrix(double elmatrix[][ELMAX1][ELMAX2], double elvector[], int anal_type, int lnflg, double time, double dt, double theta);
};
class Element24 : public Element
{
//４節点要素（４角形要素）
public:
Element24(int ID) : Element(ID)
{
node = new Node*[4];
for (int i=0; i<4; i++)
node[i] = NULL;
}
~Element24()
{
delete[] node;
}
Node* getNode(int i) { return node[i]; }
void setMaterial(Material* material) { this->material = material; }
Material* getMaterial() { return material; }
Load* getLoad() { return load; }
int getNumnode() { return 4; }
void read(string strArray[], int numat, Material** matArray, int numld, Load** loadArray, int numnp, Node** nodeArray, ifstream* fp_in);
void write(ofstream* fp_out);
//形状関数
void getRST(double rst[]);
double deriv(double r, double s, double t, double h[], double d[][3]);
//微分演算
void grad(double h[], double d[][3], double value[][3]);
//要素行列
void elMatrix(double elmatrix[][ELMAX1][ELMAX2], double elvector[], int anal_type, int lnflg, double time, double dt, double theta);
};
class Element34 : public Element
{
//４節点要素（４面体）
public:
Element34(int ID) : Element(ID)
{
node = new Node*[4];
for (int i=0; i<4; i++)
node[i] = NULL;
}
~Element34()
{
delete[] node;
}
Node* getNode(int i) { return node[i]; }
void setMaterial(Material* material) { this->material = material; }
Material* getMaterial() { return material; }
Load* getLoad() { return load; }
int getNumnode() { return 4; }
void read(string strArray[], int numat, Material** matArray, int numld, Load** loadArray, int numnp, Node** nodeArray, ifstream* fp_in);
void write(ofstream* fp_out);
//形状関数
void getRST(double rst[]);
double deriv(double r, double s, double t, double h[], double d[][3]);
//微分演算
void grad(double h[], double d[][3], double value[][3]);
//要素行列
void elMatrix(double elmatrix[][ELMAX1][ELMAX2], double elvector[], int anal_type, int lnflg, double time, double dt, double theta);
};
class Element38 : public Element
{
//８節点要素（６面体要素）
public:
Element38(int ID) : Element(ID)
{
node = new Node*[8];
for (int i=0; i<8; i++)
node[i] = NULL;
}
~Element38()
{
delete[] node;
}
Node* getNode(int i) { return node[i]; }
void setMaterial(Material* material) { this->material = material; }
Material* getMaterial() { return material; }
Load* getLoad() { return load; }
int getNumnode() { return 8; }
void read(string strArray[], int numat, Material** matArray, int numld, Load** loadArray, int numnp, Node** nodeArray, ifstream* fp_in);
void write(ofstream* fp_out);
//形状関数
void getRST(double rst[]);
double deriv(double r, double s, double t, double h[], double d[][3]);
//微分演算
void grad(double h[], double d[][3], double value[][3]);
//要素行列
void elMatrix(double elmatrix[][ELMAX1][ELMAX2], double elvector[], int anal_type, int lnflg, double time, double dt, double theta);
};
ここで定義されているメンバー関数の実装は次のようにファイル fem.cpp に与えます。
//要素クラス
//数値積分用定数
int Element::nint = 2;
double Element::xg[4][4] =
{ { 0.e0, 0.e0, 0.e0, 0.e0},
{-.57735026919e0, .57735026919e0, 0.e0, 0.e0},
{-.77459666924e0, 0.e0, .77459666924e0, 0.e0},
{-.86113631159e0,-.33998104358e0, .33998104358e0, .86113631159e0} };
double Element::wgt[4][4] =
{ { 2.e0, 0.e0, 0.e0, 0.e0},
{1.00000000000e0,1.00000000000e0, 0.e0, 0.e0},
{ .55555555556e0, .88888888889e0, .55555555556e0, 0.e0},
{ .34785484514e0, .65214515486e0, .65214515486e0, .34785484514e0} };
Element* Element::getElement(int ID, int type)
{
if(type == 23)
{
return new Element23(ID);
}
else if(type == 24)
{
return new Element24(ID);
}
else if(type == 34)
{
return new Element34(ID);
}
else if(type == 38)
{
return new Element38(ID);
}
else
{
cout << "サポートされていない要素タイプ " << type << " がありました！" << endl;
return NULL;
}
}
void Element::read(int numnode, string strArray[], int numat, Material** matArray, int numld, Load** loadArray, int numnp, Node** nodeArray, ifstream* fp_in)
{
char str[MAX_COLUMN];
int matID,loadID,nodeID[3];
//ID、物性番号、荷重番号、節点番号
strcpy(str,strArray[0].c_str()); ID = atoi(str);
strcpy(str,strArray[1].c_str()); type = atoi(str);
strcpy(str,strArray[2].c_str()); matID = atoi(str);
strcpy(str,strArray[3].c_str()); loadID = atoi(str);
for (int n=0; ngetID() == matID)
{
material = matArray[n];
find = true;
break;
}
}
if(!find)
cout << "要素番号" << ID <<" の物性番号" << matID << " がありません！" << endl;
load = NULL;
for (int n=0; n<numld; n++)
{
if(loadArray[n]->getID() == loadID)
{
load = loadArray[n];
break;
}
}
//この方法は節点数が多くなると効率が悪い
for (int i=0; i<numnode; i++)
{
find = false;
for (int n=0; n<numnp; n++)
{
if(nodeArray[n]->getID() == nodeID[i])
{
node[i] = nodeArray[n];
find = true;
break;
}
}
if(!find)
cout << "要素番号" << ID <<" の節点番号" << nodeID[i] << " がありません！" << endl;
}
}
void Element::write(int numnode, ofstream* fp_out)
{
int loadID = 0;
if(load != NULL)
loadID = load->getID();
*fp_out << ID << "," << type << "," << material->getID() << "," << loadID << ",";
for (int i=0; i<numnode; i++)
*fp_out << node[i]->getID() << ",";
*fp_out << endl;
}
void Element23::read(string strArray[], int numat, Material** matArray, int numld, Load** loadArray, int numnp, Node** nodeArray, ifstream* fp_in)
{
Element::read(3, strArray, numat, matArray, numld, loadArray, numnp, nodeArray, fp_in);
}
void Element23::write(ofstream* fp_out)
{
Element::write(3, fp_out);
}
void Element23::getRST(double rst[])
{
//deriv()をよぶときの標準座標
for (int i=0; i<3; i++)
{
rst[i] = 0.0;
for (int n=0; n<3; n++)
rst[i] += node[n]->getCoord(i)/3.0;
}
}
double Element23::deriv(double x, double y, double z, double h[], double d[][3])
{
double v1[3], v2[3], v3[3];
double v12[3],v23[3],v31[3];
double s1[3], s2[3], s3[3];
double s[3];
v1[0] = x - node[0]->getCoord(0);
v1[1] = y - node[0]->getCoord(1);
v1[2] = z - node[0]->getCoord(2);
v2[0] = x - node[1]->getCoord(0);
v2[1] = y - node[1]->getCoord(1);
v2[2] = z - node[1]->getCoord(2);
v3[0] = x - node[2]->getCoord(0);
v3[1] = y - node[2]->getCoord(1);
v3[2] = z - node[2]->getCoord(2);
for (int i=0; i<3; i++)
{
v12[i] = v1[i] - v2[i];
v23[i] = v2[i] - v3[i];
v31[i] = v3[i] - v1[i];
}
Element::outer(v2,v3,s1);
Element::outer(v3,v1,s2);
Element::outer(v1,v2,s3);
Element::outer(v12,v23,s);
double dets = sqrt(Element::inner(s,s));
double det = 0.5*dets;
double dets2 = dets*dets;
if(det > 1.0e-12)
{
h[0] = Element::inner(s1,s)/dets2;
h[1] = Element::inner(s2,s)/dets2;
h[2] = Element::inner(s3,s)/dets2;
}
else
{
cout << ID << "番目の要素をチェックしてください" << endl;
for (int i=0; i<3; i++)
h[i] = 0.0;
}
double dh1[3],dh2[3],dh3[3];
Element::outer(s,v23,dh1);
Element::outer(s,v31,dh2);
Element::outer(s,v12,dh3);
if(det > 1.0e-12)
{
for (int i=0; i<3; i++)
{
d[0][i] = dh1[i]/dets2;
d[1][i] = dh2[i]/dets2;
d[2][i] = dh3[i]/dets2;
}
}
return det;
}
void Element23::grad(double h[], double d[][3], double value[][3])
{
for (int i=0; i<3; i++)
{
value[0][i] = 0.0;
value[1][i] = 0.0;
}
for (int n=0; n<3; n++)
{
if(node[n]->getValue() != NULL)
{
double nodevalue1 = node[n]->getValue(0);
double nodevalue2 = node[n]->getValue(1);
for (int i=0; i<3; i++)
{
value[0][i] += d[n][i]*nodevalue1;
value[1][i] += d[n][i]*nodevalue2;
}
}
}
}
void Element23::elMatrix(double elmatrix[][ELMAX1][ELMAX2], double elvector[], int anal_type, int lnflg, double time, double dt, double theta)
{
}
void Element24::read(string strArray[], int numat, Material** matArray, int numld, Load** loadArray, int numnp, Node** nodeArray, ifstream* fp_in)
{
Element::read(4, strArray, numat, matArray, numld, loadArray, numnp, nodeArray, fp_in);
}
void Element24::write(ofstream* fp_out)
{
Element::write(4, fp_out);
}
void Element24::getRST(double rst[])
{
//deriv()をよぶときの標準座標
for (int i=0; i<3; i++)
rst[i] = 0.0;
}
double Element24::deriv(double r, double s, double t, double h[], double d[][3])
{
double rp,sp,rm,sm;
double dh[2][4];
double xj[3][3];
double xji[3][3];
rp = 1.0 + r; sp = 1.0 + s; rm = 1.0 - r; sm = 1.0 - s;
h[0] = 0.25*rm*sm; h[1] = 0.25*rp*sm; h[2] = 0.25*rp*sp; h[3] = 0.25*rm*sp;
dh[0][0] = -0.25*sm; dh[0][1] = 0.25*sm; dh[0][2] = 0.25*sp; dh[0][3] = -0.25*sp;
dh[1][0] = -0.25*rm; dh[1][1] = -0.25*rp; dh[1][2] = 0.25*rp; dh[1][3] = 0.25*rm;
for (int i=0; i<3; i++)
{
xj[0][i] = 0.0;
xj[1][i] = 0.0;
for (int j=0; j<4; j++)
{
xj[0][i] += dh[0][j]*node[j]->getCoord(i);
xj[1][i] += dh[1][j]*node[j]->getCoord(i);
}
}
Element::outer(xj[0],xj[1],xj[2]);
double det = sqrt(xj[2][0]*xj[2][0] + xj[2][1]*xj[2][1] + xj[2][2]*xj[2][2]);
if(det > ZERO_EPS)
{
for (int i=0; i<3; i++)
xj[2][i] = xj[2][i]/det; //xj[2]はこの面の単位法線べクトル
}
else
{
cout << ID << "番目の要素をチェックしてください" << endl;
}
xji[0][0] = (xj[1][1]*xj[2][2] - xj[1][2]*xj[2][1])/det;
xji[0][1] = (xj[0][2]*xj[2][1] - xj[0][1]*xj[2][2])/det;
xji[0][2] = (xj[0][1]*xj[1][2] - xj[0][2]*xj[1][1])/det;
xji[1][0] = (xj[1][2]*xj[2][0] - xj[1][0]*xj[2][2])/det;
xji[1][1] = (xj[0][0]*xj[2][2] - xj[2][0]*xj[0][2])/det;
xji[1][2] = (xj[1][0]*xj[0][2] - xj[1][2]*xj[0][0])/det;
xji[2][0] = (xj[2][1]*xj[1][0] - xj[2][0]*xj[1][1])/det;
xji[2][1] = (xj[2][0]*xj[0][1] - xj[2][1]*xj[0][0])/det;
xji[2][2] = (xj[0][0]*xj[1][1] - xj[0][1]*xj[1][0])/det;
for (int j=0; j<4; j++)
{
for (int i=0; i<3; i++)
{
d[j][i] = 0.0;
for (int k=0; k<2; k++)
d[j][i] += xji[i][k]*dh[k][j];
}
}
return det;
}
void Element24::grad(double h[], double d[][3], double value[][3])
{
for (int i=0; i<3; i++)
{
value[0][i] = 0.0;
value[1][i] = 0.0;
}
for (int n=0; n<4; n++)
{
if(node[n]->getValue() != NULL)
{
double nodevalue1 = node[n]->getValue(0);
double nodevalue2 = node[n]->getValue(1);
for (int i=0; i<3; i++)
{
value[0][i] += d[n][i]*nodevalue1;
value[1][i] += d[n][i]*nodevalue2;
}
}
}
}
void Element24::elMatrix(double elmatrix[][ELMAX1][ELMAX2], double elvector[], int anal_type, int lnflg, double time, double dt, double theta)
{
//４節点要素（勾配場)
double h[4];
double d[4][3];
for (int ntype=0; ntype<3; ntype++)
{
for (int i=0; i<4; i++)
{
for (int j=0; j<4; j++)
elmatrix[ntype][i][j] = 0.0;
}
}
for (int i=0; i<4; i++)
elvector[i] = 0.0;
double lvalue = 0.0;
if(load != NULL)
{
double factor = 1.0;
if(load->getTable() != NULL)
{
Table* table = load->getTable();
if(table->getType() == 2) //時間テーブル
{
double dydx;
double factor1;
if(time - dt > -ZERO_EPS)
factor1 = table->getValue(time-dt, dydx);
else
factor1 = 0.0;
double factor2 = table->getValue(time, dydx);
factor = (1.0-theta)*factor1 + theta*factor2;
}
}
else
{
//時間に依存しない場合はtime=0で factor = theta
if(time < ZERO_EPS)
factor = theta;
}
lvalue = load->getValue(0)*factor;
}
for (int lx=0; lx<nint; lx++)
{
double r = xg[nint-1][lx];
for (int ly=0; ly<nint; ly++)
{
double s = xg[nint-1][ly];
double t = 0.0;
double det = deriv(r, s, t, h, d);
double wt = wgt[nint-1][lx]*wgt[nint-1][ly];
double wtdet = wt*det;
double array[3][3];
for (int i=0; i<3; i++)
{
for (int j=0; j<3; j++)
{
if(i == j)
array[i][j] = material->getValue(0);
else
array[i][j] = 0.0;
}
}
double dvect[3][3] = {{1.0, 0.0, 0.0},{0.0, 1.0, 0.0},{0.0, 0.0, 1.0}};
//物性が非線形の場合
if(material->getTable(0) != NULL)
{
double v[3] = { 0.0, 0.0, 0.0 };
double gradv[2][3];
grad(h, d, gradv);
for (int i=0; i<3; i++)
v[i] = -((1.0-theta)*gradv[1][i] + theta*gradv[0][i]);
double absv = sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
double dydx;
double f = material->getTable(0)->getValue(absv, dydx);
if(absv > ZERO_EPS)
{
for (int i=0; i<3; i++)
{
for (int j=0; j<3; j++)
{
if(i == j)
{
array[i][j] = f;
dvect[i][j] = f + dydx*v[i]*v[j]/absv;
}
else
{
array[i][j] = 0.0;
dvect[i][j] = dydx*v[i]*v[j]/absv;
}
}
}
}
else
{
for (int i=0; i<3; i++)
{
for (int j=0; j<3; j++)
{
if(i == j)
{
array[i][j] = f;
dvect[i][j] = f;
}
else
{
array[i][j] = 0.0;
dvect[i][j] = 0.0;
}
}
}
}
}
if(anal_type == 1)
{
//軸対称
double R = 0.0;
for (int i=0; i<4; i++)
R += h[i]*node[i]->getCoord(0);
double rwtdet = 2.0*PI*R*wtdet;
for (int i=0; i<4; i++)
{
for (int j=i; j<4; j++)
{
elmatrix[0][i][j] += h[i]*material->getValue(1)*h[j]*rwtdet;
elmatrix[1][i][j] += (d[i][0]*array[0][0]*d[j][0] + d[i][0]*array[0][1]*d[j][1]
+ d[i][1]*array[1][0]*d[j][0] + d[i][1]*array[1][1]*d[j][1])*rwtdet;
elmatrix[2][i][j] += (d[i][0]*dvect[0][0]*d[j][0] + d[i][0]*dvect[0][1]*d[j][1]
+ d[i][1]*dvect[1][0]*d[j][0] + d[i][1]*dvect[1][1]*d[j][1])*rwtdet;
}
}
if(load != NULL)
{
for (int i=0; i<4; i++)
elvector[i] += h[i]*lvalue*rwtdet;
}
}
else if(anal_type == 2)
{
//２次元
for (int i=0; i<4; i++)
{
for (int j=i; j<4; j++)
{
elmatrix[0][i][j] += h[i]*material->getValue(1)*h[j]*wtdet;
elmatrix[1][i][j] += (d[i][0]*array[0][0]*d[j][0] + d[i][0]*array[0][1]*d[j][1]
+ d[i][1]*array[1][0]*d[j][0] + d[i][1]*array[1][1]*d[j][1])*wtdet;
elmatrix[2][i][j] += (d[i][0]*dvect[0][0]*d[j][0] + d[i][0]*dvect[0][1]*d[j][1]
+ d[i][1]*dvect[1][0]*d[j][0] + d[i][1]*dvect[1][1]*d[j][1])*wtdet;
}
}
if(load != NULL)
{
for (int i=0; i<4; i++)
elvector[i] += h[i]*lvalue*wtdet;
}
}
else
{
//３次元
}
}
}
//コピー
for (int ntype=0; ntype<3; ntype++)
{
for (int i=0; i<4; i++)
{
for (int j=0; j<i; j++)
elmatrix[ntype][i][j] = elmatrix[ntype][j][i];
}
}
}
void Element34::read(string strArray[], int numat, Material** matArray, int numld, Load** loadArray, int numnp, Node** nodeArray, ifstream* fp_in)
{
Element::read(4, strArray, numat, matArray, numld, loadArray, numnp, nodeArray, fp_in);
}
void Element34::write(ofstream* fp_out)
{
Element::write(3, fp_out);
}
void Element34::getRST(double rst[])
{
//deriv()をよぶときの標準座標
for (int i=0; i<3; i++)
{
rst[i] = 0.0;
for (int n=0; n<4; n++)
rst[i] += node[n]->getCoord(i)/4.0;
}
}
double Element34::deriv(double x, double y, double z, double h[], double d[][3])
{
double v12[3],v13[3],v14[3],v23[3],v24[3];
double vp1[3],vp2[3],vp3[3],vp4[3];
double s[4][3];
for (int i=0; i<3; i++)
{
v12[i] = node[1]->getCoord(i) - node[0]->getCoord(i);
v13[i] = node[2]->getCoord(i) - node[0]->getCoord(i);
v14[i] = node[3]->getCoord(i) - node[0]->getCoord(i);
v23[i] = node[2]->getCoord(i) - node[1]->getCoord(i);
v24[i] = node[3]->getCoord(i) - node[1]->getCoord(i);
}
Element::outer(v24,v23,s[0]);
Element::outer(v13,v14,s[1]);
Element::outer(v14,v12,s[2]);
Element::outer(v12,v13,s[3]);
double det = -Element::inner(s[0],v12)/6.0; //det = △/6
if(det < SMALL_EPS*SMALL_EPS*SMALL_EPS)
{
cout << ID << "番目の要素をチェックしてください" << endl;
}
vp1[0] = x - node[1]->getCoord(0);
vp1[1] = y - node[1]->getCoord(1);
vp1[2] = z - node[1]->getCoord(2);
h[0] = Element::inner(vp1,s[0])/(6.0*det);
vp2[0] = x - node[2]->getCoord(0);
vp2[1] = y - node[2]->getCoord(1);
vp2[2] = z - node[2]->getCoord(2);
h[1] = Element::inner(vp2,s[1])/(6.0*det);
vp3[0] = x - node[3]->getCoord(0);
vp3[1] = y - node[3]->getCoord(1);
vp3[2] = z - node[3]->getCoord(2);
h[2] = Element::inner(vp3,s[2])/(6.0*det);
vp4[0] = x - node[0]->getCoord(0);
vp4[1] = y - node[0]->getCoord(1);
vp4[2] = z - node[0]->getCoord(2);
h[3] = Element::inner(vp4,s[3])/(6.0*det);
for (int i=0; i<3; i++)
{
for (int j=0; j<4; j++)
d[j][i] = s[j][i]/(6.0*det);
}
return det;
}
void Element34::grad(double h[], double d[][3], double value[][3])
{
for (int i=0; i<3; i++)
{
value[0][i] = 0.0;
value[1][i] = 0.0;
}
for (int n=0; n<4; n++)
{
if(node[n]->getValue() != NULL)
{
double nodevalue1 = node[n]->getValue(0);
double nodevalue2 = node[n]->getValue(1);
for (int i=0; i<3; i++)
{
value[0][i] += d[n][i]*nodevalue1;
value[1][i] += d[n][i]*nodevalue2;
}
}
}
}
void Element34::elMatrix(double elmatrix[][ELMAX1][ELMAX2], double elvector[], int anal_type, int lnflg, double time, double dt, double theta)
{
}
void Element38::read(string strArray[], int numat, Material** matArray, int numld, Load** loadArray, int numnp, Node** nodeArray, ifstream* fp_in)
{
Element::read(8, strArray, numat, matArray, numld, loadArray, numnp, nodeArray, fp_in);
}
void Element38::write(ofstream* fp_out)
{
Element::write(8, fp_out);
}
void Element38::getRST(double rst[])
{
//deriv()をよぶときの標準座標
for (int i=0; i<3; i++)
rst[i] = 0.0;
}
double Element38::deriv(double r, double s, double t, double h[], double d[][3])
{
double rp,sp,tp,rm,sm,tm;
double dh[3][8];
double xj[3][3],xji[3][3];
rp = 1.0 + r; sp = 1.0 + s; tp = 1.0 + t;
rm = 1.0 - r; sm = 1.0 - s; tm = 1.0 - t;
h[0] = 0.125*rm*sm*tm; h[1] = 0.125*rp*sm*tm; h[2] = 0.125*rp*sp*tm; h[3] = 0.125*rm*sp*tm;
h[4] = 0.125*rm*sm*tp; h[5] = 0.125*rp*sm*tp; h[6] = 0.125*rp*sp*tp; h[7] = 0.125*rm*sp*tp;
dh[0][0] = -0.125*sm*tm; dh[1][0] = -0.125*rm*tm; dh[2][0] = -0.125*rm*sm;
dh[0][1] = -dh[0][0]; dh[1][1] = -0.125*rp*tm; dh[2][1] = -0.125*rp*sm;
dh[0][2] = 0.125*sp*tm; dh[1][2] = -dh[1][1]; dh[2][2] = -0.125*rp*sp;
dh[0][3] = -dh[0][2]; dh[1][3] = -dh[1][0]; dh[2][3] = -0.125*rm*sp;
dh[0][4] = -0.125*sm*tp; dh[1][4] = -0.125*rm*tp; dh[2][4] = -dh[2][0];
dh[0][5] = -dh[0][4]; dh[1][5] = -0.125*rp*tp; dh[2][5] = -dh[2][1];
dh[0][6] = 0.125*sp*tp; dh[1][6] = -dh[1][5]; dh[2][6] = -dh[2][2];
dh[0][7] = -dh[0][6]; dh[1][7] = -dh[1][4]; dh[2][7] = -dh[2][3];
for (int i=0; i<3; i++)
{
for (int j=0; j<3; j++)
{
xj[i][j] = 0.0;
for (int k=0; k<8; k++)
{
xj[i][j] += dh[i][k]*node[k]->getCoord(j);
}
}
}
double det = xj[0][0]*(xj[1][1]*xj[2][2] - xj[1][2]*xj[2][1])
+ xj[0][1]*(xj[1][2]*xj[2][0] - xj[1][0]*xj[2][2])
+ xj[0][2]*(xj[1][0]*xj[2][1] - xj[1][1]*xj[2][0]);
if(det < SMALL_EPS*SMALL_EPS*SMALL_EPS)
{
cout << ID << "番目の要素をチェックしてください" << endl;
}
xji[0][0] = (xj[1][1]*xj[2][2] - xj[1][2]*xj[2][1])/det;
xji[0][1] = (xj[0][2]*xj[2][1] - xj[0][1]*xj[2][2])/det;
xji[0][2] = (xj[0][1]*xj[1][2] - xj[0][2]*xj[1][1])/det;
xji[1][0] = (xj[1][2]*xj[2][0] - xj[1][0]*xj[2][2])/det;
xji[1][1] = (xj[0][0]*xj[2][2] - xj[2][0]*xj[0][2])/det;
xji[1][2] = (xj[1][0]*xj[0][2] - xj[1][2]*xj[0][0])/det;
xji[2][0] = (xj[2][1]*xj[1][0] - xj[2][0]*xj[1][1])/det;
xji[2][1] = (xj[2][0]*xj[0][1] - xj[2][1]*xj[0][0])/det;
xji[2][2] = (xj[0][0]*xj[1][1] - xj[0][1]*xj[1][0])/det;
for (int j=0; j<8; j++)
{
for (int i=0; i<3; i++)
{
d[j][i] = 0.0;
for (int k=0; k<3; k++)
d[j][i] += xji[i][k]*dh[k][j];
}
}
return det;
}
void Element38::grad(double h[], double d[][3], double value[][3])
{
for (int i=0; i<3; i++)
{
value[0][i] = 0.0;
value[1][i] = 0.0;
}
for (int n=0; n<8; n++)
{
if(node[n]->getValue() != NULL)
{
double nodevalue1 = node[n]->getValue(0);
double nodevalue2 = node[n]->getValue(1);
for (int i=0; i<3; i++)
{
value[0][i] += d[n][i]*nodevalue1;
value[1][i] += d[n][i]*nodevalue2;
}
}
}
}
void Element38::elMatrix(double elmatrix[][ELMAX1][ELMAX2], double elvector[], int anal_type, int lnflg, double time, double dt, double theta)
{
}
ここでは、２次元の四角形要素以外の実装は長くなるので一部しか行っていません。後はこれを参考にして作成してください。
また後で使うのは四角形要素だけなので記述するのが面倒な場合はそれ以外の実装部分を省略してもかまいません。