コード:
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<fstream>
#include<iostream>
#include<iomanip>
#include <time.h>
using namespace std;
/*************************************************************************
4節点四角形要素についてのプログラム
2次元定常熱伝導方程式
reference code "Visual Basicによる3次元熱伝導解析プログラム" ,黒田英夫 著
*************************************************************************/
//標準モジュールに用いる変数
int KnumMax, Mprop, Mstat, Nall, NboundP, Ndim, Nelem, Ngaus, Nmat, NMdivMax, Npoint, Nvar;
int *KnumPoint, *KnumPR, *NMdiv, **NNelem, *NPrefer, *NrestP, *NumMat;
double **Cspec, *PosGaus, **Rlamda, *Rou, *TempFixP, *TempInit, **TempMat, *Tu, *WeightG, *Xdata, *Ydata, *XE, *YE;
//フォームモジュールに用いる変数
int Nbound, NumEff, istep, mb;
double timst;
int *Jcolnum, *Krowpos, *Nrest;
double *Dd, *Ds, *Dt, *Dx, *Dy, **Ek,**RinvJ, *RJmat, *Rld, *Shape, *TE, *TempFix, *Tf, *TG, *Tk, *TUx, *Ua;
void variable(){
int maxnode = 8;
int maxnode2 = maxnode + 5;
Ngaus = 2;
Ndim = 1;
//標準モジュールの変数
//single pointer int form
NPrefer = new int[maxnode2]();
//double pointer int form
//single pointer double form
PosGaus = new double[2]();
Rou = new double[1]();
WeightG = new double[2]();
//double pointer double form
Cspec = new double*[1];
Rlamda = new double*[1];
TempMat = new double*[2];
for (int i = 1; i <= maxnode2; i++){
Cspec[i] = new double[2]();
Rlamda[i] = new double[2]();
TempMat[i] = new double[2]();
}
//フォームモジュールの変数
//single pointer int form
Jcolnum = new int[maxnode2]();
Krowpos = new int[maxnode2]();
//double pointer int form
//single pointer double form
Shape = new double[4]();
Ds = new double[4]();
Dt = new double[4]();
TE = new double[4]();
XE = new double[4]();
YE = new double[4]();
Dx = new double[4]();
Dy = new double[4]();
TUx = new double[maxnode2]();
RJmat = new double[2]();
Rld = new double[maxnode2]();
//double pointer double form
Ek = new double*[maxnode2];
RinvJ = new double*[2];
for (int i = 1; i <= 2; i++){
RinvJ[i] = new double[2]();
}
for (int i = 1; i <= maxnode2; i++){
Ek[i] = new double[maxnode2]();
}
}
void input(){
//INITIAL
ifstream fin4("initial.dat");
if (!fin4){
cout << "File for saving data cannot be opened !" << endl;
}
/********************************************************************************
initial.dat ファイルの中身
48 35 1 1 12 4
********************************************************************************/
fin4 >> KnumMax >> Nelem >> Mprop >> Nmat >> NboundP >> Nvar;
fin4.close();
//値の設定
Npoint = KnumMax;
Nall = Npoint;
NMdivMax = Nmat;
//配列の確保
int NboundPX = Nbound + 3;
int KnumMax2 = KnumMax + 3;
int Nelem2 = Nelem + 3;
Xdata = new double[KnumMax2]();
Ydata = new double[KnumMax2]();
XE = new double[KnumMax2]();
YE = new double[KnumMax2]();
NumMat = new int[Nelem2]();
NMdiv = new int[Nelem2]();
KnumPoint = new int[KnumMax2]();
Tf = new double[KnumMax2]();
TG = new double[KnumMax2]();
Tk = new double[KnumMax2]();
Ua = new double[KnumMax2]();
NrestP = new int[NboundPX]();
Nrest = new int[NboundPX]();
TempFix = new double[NboundPX]();
TempFixP = new double[NboundPX]();
Tu = new double[KnumMax2]();
TempInit = new double[Nall]();
KnumPR = new int[KnumMax2]();
NNelem = new int*[Nelem2];
for (int i = 0; i <= Nelem2; i++){
NNelem[i] = new int[4]();
}
//節点の入力
ifstream fin("node.dat");
if (!fin){
cout << "File for saving data cannot be opened !" << endl;
}
/********************************************************************************
node.dat ファイルの中身
1 0 1
2 0.1428571429 1
3 0.2857142857 1
4 0.4285714286 1
5 0.5714285714 1
6 0.7142857143 1
7 0.8571428571 1
8 1 1
9 0 0.8
10 0.1428571429 0.8
11 0.2857142857 0.8
12 0.4285714286 0.8
13 0.5714285714 0.8
14 0.7142857143 0.8
15 0.8571428571 0.8
16 1 0.8
17 0 0.6
18 0.1428571429 0.6
19 0.2857142857 0.6
20 0.4285714286 0.6
21 0.5714285714 0.6
22 0.7142857143 0.6
23 0.8571428571 0.6
24 1 0.6
25 0 0.4
26 0.1428571429 0.4
27 0.2857142857 0.4
28 0.4285714286 0.4
29 0.5714285714 0.4
30 0.7142857143 0.4
31 0.8571428571 0.4
32 1 0.4
33 0 0.2
34 0.1428571429 0.2
35 0.2857142857 0.2
36 0.4285714286 0.2
37 0.5714285714 0.2
38 0.7142857143 0.2
39 0.8571428571 0.2
40 1 0.2
41 0 0
42 0.1428571429 0
43 0.2857142857 0
44 0.4285714286 0
45 0.5714285714 0
46 0.7142857143 0
47 0.8571428571 0
48 1 0
********************************************************************************/
for (int i = 1; i <= KnumMax; i++){
fin >> KnumPoint[i] >> Xdata[i] >> Ydata[i];
}
fin.close();
//今回の場合は物質が1種類だけなので、アセンブリとかがないため局所的な節点番号と全体での節点番号は同じと見なせる。
for (int i = 1; i <= KnumMax; i++){
KnumPR[i] = i;
}
//要素の入力
ifstream fin2("element.dat");
if (!fin2){
cout << "File for saving data cannot be opened !" << endl;
}
/********************************************************************************
element.dat ファイルの中身
1 9
1 10
1 2
1 1
2 10
2 11
2 3
2 2
3 11
3 12
3 4
3 3
4 12
4 13
4 5
4 4
5 13
5 14
5 6
5 5
6 14
6 15
6 7
6 6
7 15
7 16
7 8
7 7
8 17
8 18
8 10
8 9
9 18
9 19
9 11
9 10
10 19
10 20
10 12
10 11
11 20
11 21
11 13
11 12
12 21
12 22
12 14
12 13
13 22
13 23
13 15
13 14
14 23
14 24
14 16
14 15
15 25
15 26
15 18
15 17
16 26
16 27
16 19
16 18
17 27
17 28
17 20
17 19
18 28
18 29
18 21
18 20
19 29
19 30
19 22
19 21
20 30
20 31
20 23
20 22
21 31
21 32
21 24
21 23
22 33
22 34
22 26
22 25
23 34
23 35
23 27
23 26
24 35
24 36
24 28
24 27
25 36
25 37
25 29
25 28
26 37
26 38
26 30
26 29
27 38
27 39
27 31
27 30
28 39
28 40
28 32
28 31
29 41
29 42
29 34
29 33
30 42
30 43
30 35
30 34
31 43
31 44
31 36
31 35
32 44
32 45
32 37
32 36
33 45
33 46
33 38
33 37
34 46
34 47
34 39
34 38
35 47
35 48
35 40
35 39
********************************************************************************/
for (int i = 1; i <= Nelem; i++){
int ie;
fin2 >> ie >> NNelem[i][1];
fin2 >> ie >> NNelem[i][2];
fin2 >> ie >> NNelem[i][3];
fin2 >> ie >> NNelem[i][4];
//for (int j = 1; j <= 4; j++){
// cout << ie << " " << NNelem[i][j] << endl;
//}
}
fin2.close();
//境界条件の入力
ifstream fin3("boundary.dat");
if (!fin3){
cout << "File for saving data cannot be opened !" << endl;
}
/********************************************************************************
boundary.dat ファイルの中身
1 1 100
2 9 100
3 17 100
4 25 100
5 33 100
6 41 100
7 8 0
8 16 0
9 24 0
10 32 0
11 40 0
12 48 0
********************************************************************************/
for (int i = 1; i <=NboundP; i++){
int ib;
fin3 >> ib >> NrestP[i] >>TempFixP[i];
}
fin3.close();
//熱伝導率の入力
ifstream fin5("conductivity.dat");
if (!fin5){
cout << "File for saving data cannot be opened !" << endl;
}
/********************************************************************************
conductivity.dat ファイルの中身
1 1 1 1 1 1
********************************************************************************/
int idum;
//初めから順に 材料種類 X方向熱伝導率 Y方向熱伝導率 密度 X方向比熱 Y方向比熱 とする。
fin5 >> idum >> Rlamda[1][1] >> Rlamda[1][2] >> Rou[1] >> Cspec[1][1] >> Cspec[1][2];
fin5.close();
//要素の材料番号の入力
ifstream fin6("material.dat");
if (!fin6){
cout << "File for saving data cannot be opened !" << endl;
}
/********************************************************************************
material.dat ファイルの中身
1 1 1
2 1 1
3 1 1
4 1 1
5 1 1
6 1 1
7 1 1
8 1 1
9 1 1
10 1 1
11 1 1
12 1 1
13 1 1
14 1 1
15 1 1
16 1 1
17 1 1
18 1 1
19 1 1
20 1 1
21 1 1
22 1 1
23 1 1
24 1 1
25 1 1
26 1 1
27 1 1
28 1 1
29 1 1
30 1 1
31 1 1
32 1 1
33 1 1
34 1 1
35 1 1
********************************************************************************/
for (int i = 1; i <= Nelem; i++){
int idum2;
fin6 >> idum2 >> NMdiv[i] >> NumMat[i];
}
fin6.close();
//温度上限は200℃にしておく。(適当。)
for (int j = 1; j <= Nelem; j++){
for (int i = 1; i <= NMdiv[j]; i++){
TempMat[NumMat[i]][i] = 200;
}
}
}
void ShapeFunc(double s, double t){
/***************************************************************************************************************
形状関数計算のプログラム
shape[] 形状関数
ds[] 形状関数のξの偏微分
dt[] 形状関数のηの偏微分
***************************************************************************************************************/
double s1, s1m, t1, t1m;
s1 = 1 + s;
s1m = 1 - s;
t1 = 1 + t;
t1m = 1 - t;
Shape[1] = s1m*t1m / 4.0;
Shape[2] = s1*t1m / 4.0;
Shape[3] = s1*t1 / 4.0;
Shape[4] = s1m*t1 / 4.0;
Ds[1] = -t1m / 4.0;
Ds[2] = t1m / 4.0;
Ds[3] = t1 / 4.0;
Ds[4] = -t1 / 4.0;
Dt[1] = -s1m / 4.0;
Dt[2] = -s1 / 4.0;
Dt[3] = s1 / 4.0;
Dt[4] = s1m / 4.0;
}
void JacobElem(int ie, int numg, double dj){
/***************************************************************************************************************
要素のヤコビアンマトリックス計算のプログラム
ie 要素番号
numg ガウス点番号
dj ヤコビアンマトリックスの行列式
rinvj[][] ヤコビアンマトリックスの逆行列
rjmat[][] ヤコビアンマトリックス
tg[] ガウス点の温度
***************************************************************************************************************/
int i, ii, j;
double b1;
double **RinvJ, **RJmat;
RinvJ = new double *[2];
RJmat = new double *[2];
for (i = 1; i <= 2; i++){
RinvJ[i] = new double[2]();
RJmat[i] = new double[2]();
}
for (ii = 1; ii <= 4; ii++){
TG[numg] = TG[numg] + TE[ii] * Shape[ii];
for (j = 1; j <= 2; j++){
if (j == 1){
b1 = XE[ii];
}
else if (j == 2){
b1 = YE[ii];
}
RJmat[1][j] = RJmat[i][j] + Ds[ii] * b1;
RJmat[2][j] = RJmat[2][j] + Dt[ii] * b1;
}
}
dj = RJmat[1][1] * RJmat[2][2] - RJmat[2][1] * RJmat[1][2];
RinvJ[1][1] = RJmat[2][2] / dj;
RinvJ[1][2] = -RJmat[1][2] / dj;
RinvJ[2][1] = -RJmat[2][1] / dj;
RinvJ[2][2] = RJmat[1][1] / dj;
b1 = RJmat[1][1] * RinvJ[1][1] + RJmat[1][2] * RinvJ[2][1];
if (dj <= 0.0){
cout << "Metrics error" << endl;
cout << "Element" << ie << "Gaussian Point" << numg << "Yacobian Volume is zero or minus." << endl;
}
else{
for (ii = 1; ii <= 4; ii++){
for (i = 1; i <= 2; i++){
b1 = (RinvJ[i][1] * Ds[ii] + RinvJ[i][2] * Dt[ii])/dj;
if (i == 1){
Dx[ii] = b1;
}
else if (i == 2){
Dy[ii] = b1;
}
}
}
}
}
void StateCalc(int kmat, int kdiv, double tempmesh){
/*********************************************************
フォームモジュールでのプログラム
物性状態の計算
kmat 材料番号
kdiv 物性区分番号
tempmesh 現在温度
*********************************************************/
for (int i = 1; i <= NMdiv[kmat]; i++){
kdiv = 1;
if (tempmesh <= TempMat[kmat][i]){
goto out;
}
}
out:
int dumend = 1;
}
void ElemMatrix(int ie){
/***************************************************************************************************************
要素マトリックス計算のプログラム
ie 要素番号
xe[ii],ye[ii] 要素の4節点のX,Y座標
NNelem[ie][ii] 要素の4節点の節点番号
KnumPR[ib] 各節点番号対応の節点順位
te[],tu[] 温度
nvar 要素の変数個数
ek[i][j] 熱伝導マトリックス{Ke}
em[i][j] 熱容量マトリックス{Ce}
mstat 区分変数。0で定常、1で非定常
nummat[ie] 各要素の材料番号、変数kmatにいれる
ngaus 1軸あたりのガウス点数、今回は3
ig,jg ξ、η座標に対応
rou[] 密度
cspec[] 比熱
rlamda[] 熱伝導率
ndim 1節点あたりの変数の個数
nprefer[j] 要素マトリックスの変数番号jごとに全体マトリックスの変数番号kをいれている
dx[i],dy[i] 形状関数のX,Yの偏微分
dv 積分体積
dj ヤコビアンマトリックスの行列式|J|
weight[] 積分点の重み係数
***************************************************************************************************************/
int ib, kmat, numg, dj, kdiv;
double b1, rc, rla, tempb, s, t, dv;
kdiv = 0;
dj = 0;
for (int ii = 1; ii <= 4; ii++){
int i;
ib = NNelem[ie][ii];
i = KnumPR[ib];
XE[ii] = Xdata[i];
YE[ii] = Ydata[i];
TE[ii] = Tu[i];
}
for (int i = 1; i <= Nvar; i++){
for (int j = 1; j <= i; j++){
Ek[i][j] = 0.0;
}
}
kmat = NumMat[ie];
numg = 0;
for (int ig = 1; ig <= Ngaus; ig++){
s = PosGaus[ig];
for (int jg = 1; jg <= Ngaus; jg++){
t = PosGaus[jg];
//形状関数の計算
ShapeFunc(s, t);
//要素のヤコビアンマトリックスの計算
JacobElem(ie, numg, dj);
dv = dj*WeightG[ig] * WeightG[jg];
tempb = TG[numg];
//物性状態の計算
StateCalc(kmat, kdiv, tempb);
rc = Rou[kmat] * Cspec[kmat][kdiv];
rla = Rlamda[kmat][kdiv];
for (int i = 1; i <= Nvar; i++){
for (int j = 1; j <= i; j++){
b1 = Dx[i] * Dx[j] + Dy[i] * Dy[j];
Ek[i][j] = Ek[i][j] + rla*b1*dv;
}
}
}
}
int j = 0;
for (int ii = 1; ii <= 4; ii++){
int k, i;
ib = NNelem[ie][ii];
i = KnumPR[ib];
k = Ndim*(i - 1);
for (int id = 1; id <= Ndim; id++){
j = j + 1;
k = k + 1;
NPrefer[j] = k;
}
}
}
void SurfaceCond(int kb, int kind, int *kbuf, double *gpos){
/***************************************************************************************************************
境界条件用の面設定のプログラム
kb 面番号
kind 面分類(1:面1)
kbuf[] 面節点の基本番号
gpos[] 積分点のξη座標
***************************************************************************************************************/
int i;
if (kb == 1){
kind = 1;
gpos[1] = -1;
gpos[2] = 1;
for (i = 1; i <= 4; i++){
kbuf[i] = i;
}
}
else{
cout << "2次元の解析ではありません。" << endl;
}
}
void JacobElem2(int ie, int numg, double dj, int *knode){
/***************************************************************************************************************
要素表面のヤコビアンマトリックス計算のプログラム
ie 要素番号
numg ガウス点番号
dj ヤコビアンマトリックスの行列式値
knode[] 節点の順位
***************************************************************************************************************/
int ii, j;
double b1;
double **RinvJ, **RJmat, *XG, *YG;
RinvJ = new double*[2];
RJmat = new double*[2];
XG = new double[numg]();
YG = new double[numg]();
for (int i = 1; i <= 2; i++){
RinvJ[i] = new double[2]();
RJmat[i] = new double[2]();
}
for (int ip = 1; ip <= 4; ip++){
ii = knode[ip];
XG[numg] = XG[numg] + XE[ii] * Shape[ii];
YG[numg] = YG[numg] + YE[ii] * Shape[ii];
}
for (int ip = 1; ip <= 4; ip++){
ii = knode[ip];
for (j = 1; j <= 2; j++){
if (j == 1){
b1 = XE[ii];
}
else if (j == 2){
b1 = YE[ii];
}
RJmat[1][j] = RJmat[1][j] + Ds[ii] * b1;
RJmat[2][j] = RJmat[2][j] + Dt[ii] * b1;
}
}
//逆ヤコビアンマトリックスの配列に、面積分で使用するための値を保存
dj = RJmat[1][1] * RJmat[2][2] - RJmat[1][2] * RJmat[2][1];
RinvJ[1][1] = RJmat[2][2] / dj;
RinvJ[1][2] = -RJmat[1][2] / dj;
RinvJ[2][1] = -RJmat[2][1] / dj;
RinvJ[2][2] = RJmat[1][1] / dj;
}
void PreVecCalc(int nx, double *rvec, double *crvec){
/***************************************************************************************************************
InvC・rvecの計算プログラム
nx ベクトル次数
rvec[] 対象ベクトル
crvec[] 結果ベクトル
***************************************************************************************************************/
int jb;
for (int i = 1; i <= nx; i++){
crvec[i] = rvec[i];
}
for (int i = 1; i <= nx; i++){
for (int kb = Krowpos[i]; kb <= Krowpos[i + 1] - 2; kb++){
jb = Jcolnum[kb];
crvec[i] = crvec[i] - Rld[kb] * crvec[jb];
}
crvec[i] = crvec[i] / Rld[Krowpos[i + 1] - 1];
}
for (int i = 1; i <= nx; i++){
crvec[i] = crvec[i] / Dd[i];
}
for (int i = nx; i >= 1; i--){
crvec[i] = crvec[i] / Rld[Krowpos[i + 1] - 1];
for (int kb = Krowpos[i]; kb <= Krowpos[i + 1] - 2; kb++){
jb = Jcolnum[kb];
crvec[jb] = crvec[jb] - Rld[kb] * crvec[i];
}
}
}
void StoreIndex(int i, int j, int kbx){
/*********************************************************
フォームモジュールでのプログラム
FEM行列i,j要素の保存配列位置の算出
i 行番号
j 列番号
kbx 1次元保存配列のインデックス番号
*********************************************************/
kbx = 0;
for (int kb = Krowpos[i]; kb <= Krowpos[i + 1] - 1; kb++){
if (Jcolnum[kb] == j){
kbx = kb;
goto outfor;
}
else if (Jcolnum[kb] > j){
goto outfor;
}
}
outfor:
int dumend = 1;
}
void ICCGsolver(int mb){
/***************************************************************************************************************
ICCG法のプログラム
ICCG・・・Incomplete Cholesky Conjugate Gradient
numeff 保存された変数の合計個数
tk[k] 左下三角配列の1次元保存配列
krowpos[i] 第i行の最初の保存変数に対するtk[]でのインデックス番号。
したがって、第i行の対角要素インデックスはkrowpos[i+1]-1で得られる。
また、最終をn行としてkrowpos[n+1]にその次のインデックス番号を与えておく。
jcolnum[k] tk[k]の対応する左下三角行列の列番号
mb バンドマトリックスのハーフバンド幅
dd[Nall] 対角行列D
rld[numeff] 左下三角行列L
sd[Nall] スケーリングの対角行列S
apvec[Nall] 行列計算結果Ap
crvec[Nall] 行列計算結果C^-1r
rvec[Nall] ベクトルr
pvec[Nall] ベクトルp
cnorm 右辺ベクトルbのノルム
alp α
beta β
***************************************************************************************************************/
int kb, jb, jj, kbb, kb1, j1, j2, nk, ic;
kbb = 0;
double b1, b2, fwx, fw1, fw2, fx0, eps, epsa, beta, alp, cnorm;
fw1 = 0;
fw2 = 0;
double *sd, *apvec, *crvec, *rvec, *pvec;
sd = new double[Nall]();
apvec = new double[Nall]();
crvec = new double[Nall]();
rvec = new double[Nall]();
pvec = new double[Nall]();
for (int i = 1; i <= Nall; i++){
kb = Krowpos[i + 1] - 1;
sd[i] = 1 / sqrt(Tk[kb]);
Tk[kb] = 1;
}
for (int i = 1; i <= Nall; i++){
for (int kb = Krowpos[i]; kb <= Krowpos[i + 1] - 2; kb++){
jb = Jcolnum[kb];
Tk[kb] = Tk[kb] * sd[i] * sd[jb];
}
Tf[i] = Tf[i] * sd[i];
}
//不完全コレスキー分解
for (int i = 1; i <= Nall; i++){
b1 = 0;
for (kb = Krowpos[i]; kb <= Krowpos[i + 1] - 2; kb++){
b2 = Rld[kb];
b1 = b1 + Dd[Jcolnum[kb]] * b2*b2;
}
kb = Krowpos[i + 1] - 1;
Rld[kb] = Tk[kb] - b1;
if (Rld[kb] == 0){
cout << "calculation error" << endl;
goto end;
}
else{
Dd[i] = 1 / Rld[kb];
}
jj = mb - 1;
jb = i + jj;
if (jb > Nall){
jj = Nall - i;
}
for (int j = 1; j <= jj; j++){
jb = i + j;
//保存配列位置の算出
StoreIndex(jb, i, kbb);
if (kbb > 0 && !(Tk[kbb] == 0)){
kb1 = Krowpos[jb];
b1 = 0;
for (int kb = Krowpos[i]; kb <= Krowpos[i + 1] - 2; kb++){
j1 = Jcolnum[kb];
for (int kbx = kb1; kbx <= Krowpos[jb + 1] - 2; kbx++){
j2 = Jcolnum[kbx];
if (j1 == j2){
b1 = b1 + Dd[j1] * Rld[kb] * Rld[kbx];
kb1 = kbx + 1;
goto exit;
}
else if (j2 > j1){
goto exit;
}
}
exit:
int dummy3 = 1;
}
Rld[kbb] = Tk[kbb] - b1;
}
}
}
//定常解析場合の進行状況表示
if (Mstat == 0){
fx0 = fw1;
fwx = 0.35*fw1 + 0.65*fw2;
}
//本計算
nk = Nall / 3;
eps = 1.0E-6;
cnorm = 0;
for (int i = 1; i <= Nall; i++){
rvec[i] = Tf[i];
cnorm = cnorm + Tf[i] * Tf[i];
Tf[i] = 0;
}
cnorm = sqrt(cnorm);
epsa = eps*cnorm;
//InvC・rvecの計算
PreVecCalc(Nall, rvec, pvec);
ic = 0;
for (int ik = 1; ik <= nk; ik++){
//TimCal(1, timn);
b1 = 0;
for (int i = 1; i <= Nall; i++){
b1 = b1 + pvec[i] * rvec[i];
apvec[i] = 0;
}
for (int i = 1; i <= Nall; i++){
for (int kb = Krowpos[i]; kb <= Krowpos[i + 1] - 2; kb++){
jb = Jcolnum[kb];
apvec[i] = apvec[i] + Tk[kb] * pvec[jb];
apvec[jb] = apvec[jb] + Tk[kb] * pvec[i];
}
apvec[i] = apvec[i] + Tk[Krowpos[i + 1] - 1] * pvec[i];
}
b2 = 0;
for (int i = 1; i <= Nall; i++){
b2 = b2 + pvec[i] * apvec[i];
}
alp = b1 / b2;
for (int i = 1; i <= Nall; i++){
Tf[i] = Tf[i] + alp*pvec[i];
rvec[i] = rvec[i] - alp*apvec[i];
}
//InvC・rvecの計算
PreVecCalc(Nall, rvec, pvec);
b1 = 0;
b2 = 0;
for (int i = 1; i <= Nall; i++){
b1 = b1 + crvec[i] * apvec[i];
b2 = b2 + pvec[i] * apvec[i];
}
beta = -b1 / b2;
b1 = 0;
for (int i = 1; i <= Nall; i++){
pvec[i] = crvec[i] + beta*pvec[i];
b1 = b1 + rvec[i] * rvec[i];
}
b1 = sqrt(b1);
if (b1 <= epsa){
ic = 1;
goto exit2;
}
}
exit2:
if (ic == 0){
cout << "Solver isn't converged" << endl;
}
for (int i = 1; i <= Nall; i++){
Tf[i] = Tf[i] * sd[i];
}
end:
int dum1 = 1;
}
void KnumComp(){
/*********************************************************
標準モジュールでのプログラム
節点番号から節点順位への照合計算
*********************************************************/
int k, KnumMax;
KnumMax = 1;
for (int i = 1; i <= Npoint; i++){
k = KnumPoint[i];
if (k > KnumMax){
KnumMax = k;
}
}
double *KnumPR;
KnumPR = new double[KnumMax]();
for (int i = 1; i <= Npoint; i++){
k = KnumPoint[i];
KnumPR[k] = i;
}
}
void BandWidth(int mb){
/*********************************************************
フォームモジュールでのプログラム
バンド幅の計算
*********************************************************/
int mx,ib;
mx = 0;
ib = 0;
for (int ie = 1; ie <= Nelem; ie++){
for (int ii = 1; ii <= 3; ii++){
int j,i,k;
ib = NNelem[ie][ii];
cout << ib <<endl;
i = KnumPR[ib];
for (int jj = ii + 1; jj <= 4; jj++){
ib = NNelem[ie][jj];
j = KnumPR[ib];
k = abs(i - j);
if (k > mx){
mx = k;
}
}
}
}
mb = Ndim*(mx + 1);
}
void SurfaceSet(int ib, int kb, int *kbuf){
/*********************************************************
フォームモジュールでのプログラム
面の設定
ib 要素番号
kbuf[] 面節点の順位番号
*********************************************************/
kbuf[1] = NNelem[ib][1];
kbuf[2] = NNelem[ib][2];
kbuf[3] = NNelem[ib][3];
kbuf[4] = NNelem[ib][4];
for (int j = 1; j <= 4; j++){
kbuf[j] = KnumPR[kbuf[j]];
}
}
void InitialSet1(){
/*********************************************************
フォームモジュールでのプログラム
初期設定
*********************************************************/
Nvar = 4;
Nall = Ndim*Npoint;
//温度固定データの処理
Nbound = NboundP;
if (Nbound > 0){
int k;
k = 0;
if (NboundP > 0){
for (int i = 1; i <= NboundP; i++){
k = k + 1;
Nrest[k] = NrestP[i];
TempFix[k] = TempFixP[i];
}
}
}
for (int i = 1; i <= Nall; i++){
Tu[i] = TempInit[i];
}
PosGaus[1] = -0.57735026918962576451;
PosGaus[2] = 0.57735026918962576451;
WeightG[1] = 1.0;
WeightG[2] = 1.0;
}
void InitialSet2(){
/*********************************************************
フォームモジュールでのプログラム
初期設定2
*********************************************************/
//関連節点数の計算
int numall, nb, num,ic,ign;
int *nbuf;
nbuf = new int[8]();
numall = 0;
nb = int(Ndim*(Ndim - 1) / 2);
for (int ip = 1; ip <= Npoint; ip++){
for (int j = 1; j <= ip; j++){
nbuf[j] = 0;
}
for (int ie = 1; ie <= Nelem; ie++){
ic = 0;
for (int ii = 1; ii <= 4; ii++){
if (KnumPR[NNelem[ie][ii]] == ip){
ic = 1;
goto out;
}
}
out:
int dummy = 1;
}
num = 0;
for (int j = 1; j <= ip; j++){
if (nbuf[j] > 0){
num = num + 1;
}
}
numall = numall + num*Ndim*Ndim - nb;
}
NumEff = numall;
//保存配列との関連付け
numall = 0;
nb = 0;
for (int ip = 1; ip <= Npoint; ip++){
for (int j = 1; j <= ip; j++){
nbuf[j] = 0;
}
for (int ie = 1; ie <= Nelem; ie++){
int ic;
ic = 0;
for (int ii = 1; ii <= 4; ii++){
if (KnumPR[NNelem[ie][ii]] == ip){
ic = 1;
goto out2;
}
}
out2:
int dummy2 = 1;
if (ic > 0){
for (int ii = 1; ii <= 4; ii++){
int ib;
ib = KnumPR[NNelem[ie][ii]];
nbuf[ib] = 1;
}
}
}
for (int id = 1; id <= Ndim; id++){
int jb;
jb = 0;
nb = nb + 1;
Krowpos[nb] = numall + 1;
for (int j = 1; j <= ip; j++){
if (nbuf[j] > 0){
int k2;
if (j == ip){
k2 = id;
}
else{
k2 = Ndim;
}
for (int k = 1; k <= k2; k++){
jb = jb + 1;
numall = numall + 1;
Jcolnum[numall] = jb;
}
}
else{
jb = jb + Ndim;
}
}
}
if (ip == Npoint){
Krowpos[nb + 1] = numall + 1;
}
}
}
void RestCond(int k, int mb, double ub){
/*********************************************************
フォームモジュールでのプログラム
拘束条件処理
k 変数番号
mb ハーフバンド幅
ub 拘束値
*********************************************************/
int mx, it;
if (k > 1){
for (int kb = Krowpos[k]; kb <= Krowpos[k + 1] - 2; kb++){
if (kb > 0){
if (!ub == 0){
it = Jcolnum[kb];
Tf[it] = Tf[it] - Tk[kb] * ub;
}
Tk[kb] = 0;
}
}
}
if (k < Nall){
mx = mb;
if (k + mb - 1 > Nall){
mx = Nall - k + 1;
}
for (int j = 2; j <= mx; j++){
int kb;
kb = 0;
StoreIndex(k + j - 1, k, kb);
//保存配列位置の算出
if (kb > 0){
if (!ub == 0){
it = k + j - 1;
Tf[it] = Tf[it] - Tk[kb] * ub;
}
Tk[kb] = 0;
}
}
}
int kb;
kb = Krowpos[k + 1] - 1;
Tk[kb] = 1;
Tf[k] = ub;
}
void Store(){
/*********************************************************
フォームモジュールでのプログラム
全体マトリックスへの蓄積
*********************************************************/
for (int i = 1; i <= Nvar; i++){
for (int j = 1; i <= i; j++){
int it, jt, ib;
it = NPrefer[i];
jt = NPrefer[j];
if (it < jt){
ib = it;
it = jt;
jt = ib;
}
int kb;
kb = 0;
StoreIndex(it, jt, kb);
//保存配列位置の算出
Tk[kb] = Tk[kb] + Ek[i][j];
}
}
}
void TimCal(int Ktime, double timn$){
/*********************************************************
フォームモジュールでのプログラム
タイマー計算のためのプログラム
Ktime 時間継続指定(0:時間クリアする 1:時間継続する)
timn$ 時分秒に換算した経過時間
*********************************************************/
double rtime, rmin, rhour, sec;
rtime = 0;
//時間計測?
rtime = clock();
rtime = rtime / CLOCKS_PER_SEC;
timn$ = rtime;
if (Ktime > 0){
sec = rtime - timst;
do{
if (sec < 0){
sec = sec + 86400;
}
} while (sec < 0);
rhour = int(sec / 3600);
sec = sec - rhour * 3600;
rmin = int(sec / 60);
sec = int(sec - rmin * 60);
cout << rhour << "時間" << rmin << "分" << sec << "秒" << endl;
}
else{
timst = rtime;
}
}
void ResFile1(){
/*********************************************************
フォームモジュールでのプログラム
結果ファイルへの出力1
*********************************************************/
ofstream fout("result.dat");
if (!fout){
cout << "File for saving data cannot be opened !" << endl;
}
fout << "<<2次元定常熱伝導解析の計算結果>>" << endl;
if (Mstat == 0){
fout << "解析種類 (定常)" << endl;
}
else if (Mstat == 1){
fout << "解析種類 (非定常)" << endl;
}
fout << "------------------------------------------------------------------------" << endl;
}
void ResFile2(int mb){
/*********************************************************
フォームモジュールでのプログラム
結果ファイルへの出力2
mb ハーフバンド幅
*********************************************************/
ofstream fout("result.dat");
if (!fout){
cout << "File for saving data cannot be opened !" << endl;
}
fout << "***計算結果***" << endl;
fout << "バンド幅" << mb << endl;
fout << "保存配列のインデックス数" << NumEff << endl;
}
void ResFile3(int itime, int icycle, int istep){
/*********************************************************
フォームモジュールでのプログラム
結果ファイルへの出力3
itime 時間メッシュ番号
*********************************************************/
ofstream fout("result.dat");
if (!fout){
cout << "File for saving data cannot be opened !" << endl;
}
fout << "<節点番号,温度>" << endl;
for (int ip = 1; ip <= Npoint; ip++){
double b1, b2;
int i1, z1;
b1 = Tu[ip];
i1 = int(b1);
b2 = b1 - i1;
z1 = int(b2 * 100) / 100 + i1;
fout << KnumPoint[ip] << setw(4) << "" << z1 << endl;
}
}
void FormActivate(){
double timn$, rmesh;
int icycle, itime;
rmesh = 0;
timn$ = 0;
istep = 0;
icycle = 0;
//TimCal(0, timn$);
//cout << "計算実行中" << endl;
//cout << "しばらくお待ちください" << endl;
//cout << "" << endl;
//節点番号から節点順位への照合計算
//KnumComp();
//結果ファイルへの出力
//ResFile1();
//cout << "初期設定処理" << endl;
//TimCal(1, timn$);
//cout << "経過時間" << timn$ << endl;
//初期設定
InitialSet1();
//バンド幅の計算
BandWidth(mb);
//初期設定2
InitialSet2();
//結果ファイルへの出力2
ResFile2(mb);
itime = 0;
istep = 1;
cout << "要素計算中" << endl;
TimCal(1, timn$);
for (int ie = 1; ie <= Nelem; ie++){
//要素マトリックスの計算
ElemMatrix(ie);
//全体マトリックスへの蓄積
Store();
}
int icheck;
icheck = 0;
if (Nbound > 0){
for (int ib = 1; ib <= Nbound; ib++){
int ix, ii, k;
double ub;
ix = Nrest[ib];
ii = KnumPR[ix];
k = Ndim*(ii - 1) + 1;
ub = TempFix[ib];
//拘束条件処理
RestCond(k, mb, ub);
}
}
//ICCG法ソルバー計算
ICCGsolver(mb);
for (int i = 1; i <= Nall; i++){
Tu[i] = Tf[i];
}
//定常の場合
if (Mstat == 0){
//結果ファイルへの出力3
ResFile3(itime, icycle, istep);
//進行表示
icheck = 1;
}
TimCal(1, timn$);
cout << "計算終了" << endl;
cout << "計算結果をファイルへ出力しました。" << endl;
}
void output(){
}
int main(){
//変数の設定
variable();
//データの入力
input();
//計算の制御
FormActivate();
//計算結果の出力
output();
}