using System; using System.Collections.Generic; using System.ComponentModel; using System.Data; using System.Drawing; using System.Linq; using System.Text; using System.Windows.Forms; namespace vcsFRAME { public partial class Form1 : Form { public Form1() { InitializeComponent(); } //-------------------------------------------------------------------------------- private void Form1_Load(object sender, EventArgs e) { label1.Dock = DockStyle.Fill; label1.TextAlign = ContentAlignment.MiddleCenter; label3.Dock = DockStyle.Fill; label3.TextAlign = ContentAlignment.MiddleCenter; label5.Dock = DockStyle.Fill; label5.TextAlign = ContentAlignment.MiddleCenter; label7.Dock = DockStyle.Fill; label7.TextAlign = ContentAlignment.MiddleCenter; label9.Dock = DockStyle.Fill; label9.TextAlign = ContentAlignment.MiddleCenter; label2.Dock = DockStyle.Fill; label2.TextAlign = ContentAlignment.MiddleLeft; label4.Dock = DockStyle.Fill; label4.TextAlign = ContentAlignment.MiddleLeft; label6.Dock = DockStyle.Fill; label6.TextAlign = ContentAlignment.MiddleLeft; label8.Dock = DockStyle.Fill; label8.TextAlign = ContentAlignment.MiddleLeft; label10.Dock = DockStyle.Fill; label10.TextAlign = ContentAlignment.MiddleLeft; label1.Text = "入力ファイル"; label3.Text = "出力ファイル"; label5.Text = "開始時刻"; label7.Text = "完了時刻"; label9.Text = "計算時間"; label2.Text = ""; label4.Text = ""; label6.Text = ""; label8.Text = ""; label10.Text = ""; label11.Text = "nt,mm,ib"; } //------------------------------------------------------------------------------ private void toolStripButton1_Click(object sender, EventArgs e) { main_part(); } //-------------------------------------------------------------------------------- private void main_part() { //平面骨組構造解析 int i, j, ne, n, ia, ja; System.IO.StreamReader sr; System.IO.StreamWriter sw; string dat = ""; string[] sbuf; char delim = ','; string fnameR = ""; string fnameW = ""; string strcom = ""; int NODT; //節点総数 int NELT; //要素総数 int MATEL; //材料種別数 int KOX; //x方向変位既知節点数 int KOY; //y方向変位既知節点数 int KOZ; //回転変位既知節点数 int NF; //載荷節点数 int nt; //全自由度 int mm; //変位既知節点処理語の自由度 int ib; //バンド幅 double fact0 = 1e-30; //0判定 double ggg = 9.8; //重力加速度 //重力加速度は9.8に固定．慣性力ベクトルを求める中でgで除してgを乗じるため値は任意でかまわない DateTime sttime; //計算開始時刻 DateTime entime; //計算終了時刻 TimeSpan tspan; //計算時間 //----------------------------------------------------- //入出力データファイル指定 //----------------------------------------------------- openFileDialog1.InitialDirectory = System.IO.Directory.GetCurrentDirectory(); if (openFileDialog1.ShowDialog() == DialogResult.OK) { fnameR = openFileDialog1.FileName; } if (saveFileDialog1.ShowDialog() == DialogResult.OK) { fnameW = saveFileDialog1.FileName; } label2.Text = fnameR; label4.Text = fnameW; //----------------------------------------------------- //時間計測開始 //----------------------------------------------------- sttime = DateTime.Now; label6.Text = sttime.ToString(); //***************************************************** //データ入力 //***************************************************** sr = new System.IO.StreamReader(fnameR, System.Text.Encoding.Default); //-------------------------------------------------------------------------------------- //書出用コメント入力 //-------------------------------------------------------------------------------------- dat = sr.ReadLine(); sbuf = dat.Split(delim); strcom = sbuf[0]; //-------------------------------------------------------------------------------------- //節点数，要素数，材料種別数,X方向変位既知節点数，Y方向変位既知節点数，回転変位既知節点数，載荷点数 //-------------------------------------------------------------------------------------- dat = sr.ReadLine(); sbuf = dat.Split(delim); NODT = int.Parse(sbuf[0]); NELT = int.Parse(sbuf[1]); MATEL = int.Parse(sbuf[2]); KOX = int.Parse(sbuf[3]); KOY = int.Parse(sbuf[4]); KOZ = int.Parse(sbuf[5]); NF = int.Parse(sbuf[6]); nt = NODT * 3; //-------------------------------------------------------------------------------------- //配列寸法宣言 //-------------------------------------------------------------------------------------- int[,] kakom = new int[NELT + 1, 3]; //要素構成節点番号 double[] Em = new double[NELT + 1]; //要素弾性係数 double[] AA = new double[NELT + 1]; //要素断面積 double[] AI = new double[NELT + 1]; //要素断面二次モーメント double[] gamma = new double[NELT + 1]; //要素密度 double[] gkh = new double[NELT + 1]; //水平加速度 double[] gkv = new double[NELT + 1]; //鉛直加速度 double[] alpha = new double[NELT + 1]; //線膨張係数 int[] matno = new int[NELT + 1]; //材料種別No double[,] wm = new double[MATEL + 1, 8]; //材料種別作業用 double[] x = new double[NODT + 1]; //節点x座標 double[] y = new double[NODT + 1]; //節点y座標 double[] deltaT = new double[NODT + 1]; //節点温度変化量 int[] nokx = new int[KOX + 1]; //x方向変位既知節点番号 int[] noky = new int[KOY + 1]; //y方向変位既知節点番号 int[] nokz = new int[KOZ + 1]; //回転変位既知節点番号 int[] index = new int[nt + 1]; //変位既知節点処理のためのインデックス double[,] tsm = new double[nt + 1, nt + 1]; //全体剛性マトリックス double[] f = new double[nt + 1]; //全体節点外力ベクトル（入力値） double[] fmass = new double[nt + 1]; //全体慣性力ベクトル double[] ftemp = new double[nt + 1]; //全体温度荷重相当ベクトル double[] ftvec = new double[nt + 1]; //全体節点力ベクトル double[] fdis = new double[nt + 1]; //既知節点変位による荷重成分格納 double[] rdis = new double[nt + 1]; //既知節点変位量 double[] dis = new double[nt + 1]; //全体節点変位 double[] reac = new double[nt + 1]; //要素節点反力 double[,] sm = new double[7, 7]; //要素剛性マトリックス double[,] hen = new double[7, 7]; //座標変換行列 double[] fevec = new double[7]; //要素節点力ベクトル //------------------------------------------------------------------------------------------ //材料特性（Em,AA,AI） //------------------------------------------------------------------------------------------ for (i = 1; i <= MATEL; i++) { dat = sr.ReadLine(); sbuf = dat.Split(delim); for (j = 1; j <= 7; j++) { wm[i, j] = double.Parse(sbuf[j - 1]); } } //------------------------------------------------------------------------------------------ //要素構成節点と材料種別（node-1,node-2,matno） //------------------------------------------------------------------------------------------ for (ne = 1; ne <= NELT; ne++) { dat = sr.ReadLine(); sbuf = dat.Split(delim); kakom[ne, 1] = int.Parse(sbuf[0]); kakom[ne, 2] = int.Parse(sbuf[1]); matno[ne] = int.Parse(sbuf[2]); for (i = 1; i <= MATEL; i++) { if (i == matno[ne]) { Em[ne] = wm[i, 1]; AA[ne] = wm[i, 2]; AI[ne] = wm[i, 3]; gamma[ne] = wm[i, 4]; gkh[ne] = wm[i, 5]; gkv[ne] = wm[i, 6]; alpha[ne] = wm[i, 7]; } } } //-------------------------------------------------------------------------------------- //節点座標（x,y） //-------------------------------------------------------------------------------------- for (i = 1; i <= NODT; i++) { dat = sr.ReadLine(); sbuf = dat.Split(delim); x[i] = double.Parse(sbuf[0]); y[i] = double.Parse(sbuf[1]); deltaT[i] = double.Parse(sbuf[2]); } //-------------------------------------------------------------------------------------- //変位既知節点番号，既知節点変位量 //-------------------------------------------------------------------------------------- for (i = 1; i <= nt; i++) { rdis[i] = 0.0; } if (0 < KOX) { for (i = 1; i <= KOX; i++) { dat = sr.ReadLine(); sbuf = dat.Split(delim); nokx[i] = int.Parse(sbuf[0]); rdis[3 * nokx[i] - 2] = double.Parse(sbuf[1]); } } if (0 < KOY) { for (i = 1; i <= KOY; i++) { dat = sr.ReadLine(); sbuf = dat.Split(delim); noky[i] = int.Parse(sbuf[0]); rdis[3 * noky[i] - 1] = double.Parse(sbuf[1]); } } if (0 < KOZ) { for (i = 1; i <= KOZ; i++) { dat = sr.ReadLine(); sbuf = dat.Split(delim); nokz[i] = int.Parse(sbuf[0]); rdis[3 * nokz[i]] = double.Parse(sbuf[1]); } } //-------------------------------------------------------------------------------------- //節点番号,x方向荷重,y方向荷重,モーメント荷重 //-------------------------------------------------------------------------------------- for (i = 1; i <= nt; i++) { f[i] = 0.0; } if (0 < NF) { for (i = 1; i <= NF; i++) { dat = sr.ReadLine(); sbuf = dat.Split(delim); n = int.Parse(sbuf[0]); f[3 * n - 2] = double.Parse(sbuf[1]); f[3 * n - 1] = double.Parse(sbuf[2]); f[3 * n] = double.Parse(sbuf[3]); } } sr.Close(); //***************************************************** //入力確認出力 //***************************************************** INPDAT(fnameW, strcom, NODT, NELT, MATEL, KOX, KOY, KOZ, NF, x, y, deltaT, f, rdis, nokx, noky, nokz, kakom, Em, AA, AI, gamma, gkh, gkv, alpha, matno); //***************************************************** //メイン処理計算 //***************************************************** //----------------------------------------------------- //全体剛性行列・全体荷重ベクトルクリア //----------------------------------------------------- for (i = 1; i <= nt; i++) { ftvec[i] = 0.0; dis[i] = 0.0; reac[i] = 0.0; for (j = 1; j <= nt; j++) { tsm[i, j] = 0.0; } } //----------------------------------------------------- //要素剛性行列組立 //----------------------------------------------------- for (ne = 1; ne <= NELT; ne++) { //剛性行列組立 SM_FRM(ne, kakom, x, y, Em, AA, AI, sm); //要素剛性行列 GLTMX(ne, kakom, x, y, sm, hen); //座標変換行列と[ak]の座標変換 ASMAT(ne, 2, 3, kakom, sm, tsm); //全体行列組立：[ak]→[TAK] //慣性力ベクトル組立 MM_FRM(ne, kakom, x, y, ggg, AA, gamma, sm); //要素質量行列 GLTMX(ne, kakom, x, y, sm, hen); //座標変換行列と[am]の座標変換 MV_FRM(ne, sm, ggg, gkh, gkv, fevec); ASVEC(ne, 2, 3, kakom, fevec, fmass); //温度荷重相当ベクトル組立 TV_FRM(ne, kakom, deltaT, Em, AA, alpha, fevec, hen); ASVEC(ne, 2, 3, kakom, fevec, ftemp); } //----------------------------------------------------- //既知節点変位による荷重成分ベクトル作成 //----------------------------------------------------- for (i = 1; i <= nt; i++) { fdis[i] = 0.0; for (j = 1; j <= nt; j++) { fdis[i] = fdis[i] + rdis[j] * tsm[i, j]; } } //----------------------------------------------------- //荷重ベクトルセット //----------------------------------------------------- for (i = 1; i <= nt; i++) { ftvec[i] = f[i] + fmass[i] + ftemp[i] - fdis[i]; } //----------------------------------------------------- //変位既知節点自由度記憶 //----------------------------------------------------- for (i = 1; i <= nt; i++) { index[i] = i; reac[i] = ftvec[i]; } if (0 < KOX) { for (i = 1; i <= KOX; i++) { n = 3 * nokx[i] - 2; index[n] = 0; } } if (0 < KOY) { for (i = 1; i <= KOY; i++) { n = 3 * noky[i] - 1; index[n] = 0; } } if (0 < KOZ) { for (i = 1; i <= KOZ; i++) { n = 3 * nokz[i]; index[n] = 0; } } mm = 0; for (i = 1; i <= nt; i++) { if (index[i] > 0) { mm = mm + 1; index[mm] = i; } } //----------------------------------------------------- //変位既知節点処理（連立一次方程式縮小） //----------------------------------------------------- for (i = 1; i <= mm; i++) { ia = index[i]; reac[i] = reac[ia]; for (j = 1; j <= mm; j++) { ja = index[j]; tsm[i, j] = tsm[ia, ja]; } } //----------------------------------------------------- //バンド幅計算とフルマトリックスのバンド化記憶 //----------------------------------------------------- ib = IBAND(mm, tsm, fact0); BAND(mm, ib, tsm); //----------------------------------------------------- //対角項確認と異常終了通知 //----------------------------------------------------- for (i = 1; i <= mm; i++) { if (Math.Abs(tsm[i, 1]) < fact0) IERROR2(fnameW); } //----------------------------------------------------- //連立一次方程式（コレスキー法） //----------------------------------------------------- CHBAND10(mm, ib, tsm); //三角行列 CHBAND11(mm, ib, tsm, reac); //前進消去・後退代入 //----------------------------------------------------- //変位計算 //----------------------------------------------------- for (i = 1; i <= nt; i++) { dis[i] = 0.0; } for (i = 1; i <= mm; i++) { ia = index[i]; dis[ia] = reac[i]; //連立方程式の解としての変位セット } for (i = 1; i <= nt; i++) { reac[i] = 0.0; dis[i] = dis[i] + rdis[i]; //既知節点変位の足し込み } //----------------------------------------------------- //断面力計算 //----------------------------------------------------- for (ne = 1; ne <= NELT; ne++) { SM_FRM(ne, kakom, x, y, Em, AA, AI, sm); //要素剛性行列 GLTMX(ne, kakom, x, y, sm, hen); //座標変換行列と[ak]の座標変換 CALST_FRM(ne, kakom, dis, sm, hen, reac, tsm, Em, AA, alpha, deltaT); } //----------------------------------------------------- //時間計測完了 //----------------------------------------------------- entime = DateTime.Now; tspan = entime - sttime; label8.Text = entime.ToString(); label10.Text = "Calculation Time=" + tspan.TotalSeconds.ToString("0.00sec"); label11.Text = "nt=" + nt.ToString() + " " + "mm=" + mm.ToString() + " " + "ib=" + ib.ToString(); //***************************************************** //結果出力 //***************************************************** OUTRESUL(fnameW, NODT, NELT, reac, ftvec, x, y, dis, tsm, matno); sw = new System.IO.StreamWriter(fnameW, true, System.Text.Encoding.Default); dat = "NODT=" + NODT.ToString() + " nt=" + nt.ToString() + " mm=" + mm.ToString() + " ib=" + ib.ToString(); sw.WriteLine(dat); dat = "Calculation time=" + tspan.TotalSeconds.ToString("0.00sec"); sw.WriteLine(dat); dat = entime.ToString(); sw.WriteLine(dat); sw.Close(); MessageBox.Show("計算完了", "通知"); } //-------------------------------------------------------------------------------- private void INPDAT(string fnameW, string strcom, int NODT, int NELT, int MATEL, int KOX, int KOY, int KOZ, int NF, double[] x, double[] y, double[] deltaT, double[] f, double[] rdis, int[] nokx, int[] noky, int[] nokz, int[,] kakom, double[] Em, double[] AA, double[] AI, double[] gamma, double[] gkh, double[] gkv, double[] alpha, int[] matno) { //入力確認出力 System.IO.StreamWriter sw; string dat; int i, j, k, l, m, ne; sw = new System.IO.StreamWriter(fnameW, false, System.Text.Encoding.Default); dat = strcom; sw.WriteLine(dat); dat = "NODT,NELT,MATEL,KOX,KOY,KOZ,NF"; sw.WriteLine(dat); dat = NODT.ToString("0"); dat = dat + "," + NELT.ToString("0"); dat = dat + "," + MATEL.ToString("0"); dat = dat + "," + KOX.ToString("0"); dat = dat + "," + KOY.ToString("0"); dat = dat + "," + KOZ.ToString("0"); dat = dat + "," + NF.ToString("0"); sw.WriteLine(dat); dat = "*node characteristics"; sw.WriteLine(dat); dat = "node,x,y,fx,fy,fz,fix-x,fix-y,fix-z,rdis-x,rdis-y,rdis-z,deltaT"; sw.WriteLine(dat); for (i = 1; i <= NODT; i++) { k = 0; l = 0; m = 0; if (0 < KOX) { for (j = 1; j <= KOX; j++) { if (i == nokx[j]) k = 1; } } if (0 < KOY) { for (j = 1; j <= KOY; j++) { if (i == noky[j]) l = 1; } } if (0 < KOZ) { for (j = 1; j <= KOZ; j++) { if (i == nokz[j]) m = 1; } } dat = i.ToString("0"); dat = dat + "," + x[i].ToString("E"); dat = dat + "," + y[i].ToString("E"); dat = dat + "," + f[3 * i - 2].ToString("E"); dat = dat + "," + f[3 * i - 1].ToString("E"); dat = dat + "," + f[3 * i].ToString("E"); dat = dat + "," + k.ToString("0"); dat = dat + "," + l.ToString("0"); dat = dat + "," + m.ToString("0"); dat = dat + "," + rdis[3 * i - 2].ToString("E"); dat = dat + "," + rdis[3 * i - 1].ToString("E"); dat = dat + "," + rdis[3 * i].ToString("E"); dat = dat + "," + deltaT[i].ToString("E"); sw.WriteLine(dat); } dat = "*element characteristics"; sw.WriteLine(dat); dat = "element,node-1,node-2,E,A,I,gamma,gkh,gkv,alpha,matno"; sw.WriteLine(dat); for (ne = 1; ne <= NELT; ne++) { dat = ne.ToString("0"); dat = dat + "," + kakom[ne, 1].ToString("0"); dat = dat + "," + kakom[ne, 2].ToString("0"); dat = dat + "," + Em[ne].ToString("E"); dat = dat + "," + AA[ne].ToString("E"); dat = dat + "," + AI[ne].ToString("E"); dat = dat + "," + gamma[ne].ToString("E"); dat = dat + "," + gkh[ne].ToString("E"); dat = dat + "," + gkv[ne].ToString("E"); dat = dat + "," + alpha[ne].ToString("E"); dat = dat + "," + matno[ne].ToString("0"); sw.WriteLine(dat); } sw.Close(); } //-------------------------------------------------------------------------------- private void OUTRESUL(string fnameW, int NODT, int NELT, double[] reac, double[] ftvec, double[] x, double[] y, double[] dis, double[,] tsm, int[] matno) { //結果出力（Appendモード） int i, ne; System.IO.StreamWriter sw; string dat; sw = new System.IO.StreamWriter(fnameW, true, System.Text.Encoding.Default); dat = "*displacements and forces"; sw.WriteLine(dat); dat = "node,x-coord,y-coord,x-dir,y-dir,z-dir,x-reac,y-reac,z-reac,ftvec-x,ftvec-y,ftvec-z"; sw.WriteLine(dat); for (i = 1; i <= NODT; i++) { dat = i.ToString("0"); dat = dat + "," + x[i].ToString("E"); dat = dat + "," + y[i].ToString("E"); dat = dat + "," + dis[3 * i - 2].ToString("E"); dat = dat + "," + dis[3 * i - 1].ToString("E"); dat = dat + "," + dis[3 * i].ToString("E"); dat = dat + "," + reac[3 * i - 2].ToString("E"); dat = dat + "," + reac[3 * i - 1].ToString("E"); dat = dat + "," + reac[3 * i].ToString("E"); dat = dat + "," + ftvec[3 * i - 2].ToString("E"); dat = dat + "," + ftvec[3 * i - 1].ToString("E"); dat = dat + "," + ftvec[3 * i].ToString("E"); sw.WriteLine(dat); } dat = "*stress resultants"; sw.WriteLine(dat); dat = "element,Ni,Si,Mi,Nj,Sj,Mj"; sw.WriteLine(dat); for (ne = 1; ne <= NELT; ne++) { dat = ne.ToString("0"); for (i = 1; i <= 6; i++) { dat = dat + "," + tsm[ne, i].ToString("E"); } sw.WriteLine(dat); } sw.Close(); } //-------------------------------------------------------------------------------- private void SM_FRM(int ne, int[,] kakom, double[] x, double[] y, double[] Em, double[] AA, double[] AI, double[,] ak) { //要素剛性行列 int i, j; double AL, EA, EI; //行列初期化 for (i = 1; i <= 6; i++) { for (j = 1; j <= 6; j++) { ak[i, j] = 0.0; } } //要素長さ i = kakom[ne, 1]; j = kakom[ne, 2]; AL = Math.Sqrt((x[j] - x[i]) * (x[j] - x[i]) + (y[j] - y[i]) * (y[j] - y[i])); //要素剛性行列 EA = Em[ne] * AA[ne]; EI = Em[ne] * AI[ne]; ak[1, 1] = EA / AL; ak[2, 1] = 0.0; ak[3, 1] = 0.0; ak[4, 1] = -EA / AL; ak[5, 1] = 0.0; ak[6, 1] = 0.0; ak[2, 2] = 12.0 * EI / AL / AL / AL; ak[3, 2] = 6.0 * EI / AL / AL; ak[4, 2] = 0.0; ak[5, 2] = -12.0 * EI / AL / AL / AL; ak[6, 2] = 6.0 * EI / AL / AL; ak[3, 3] = 4.0 * EI / AL; ak[4, 3] = 0.0; ak[5, 3] = -6.0 * EI / AL / AL; ak[6, 3] = 2.0 * EI / AL; ak[4, 4] = EA / AL; ak[5, 4] = 0.0; ak[6, 4] = 0.0; ak[5, 5] = 12.0 * EI / AL / AL / AL; ak[6, 5] = -6.0 * EI / AL / AL; ak[6, 6] = 4.0 * EI / AL; //対称項セット for (i = 1; i <= 5; i++) { for (j = i + 1; j <= 6; j++) { ak[i, j] = ak[j, i]; } } } //------------------------------------------------------------------------------ private void MM_FRM(int ne, int[,] kakom, double[] x, double[] y, double ggg, double[] AA, double[] gamma, double[,] am) { //要素質量行列 int i, j; double AL; //行列初期化 for (i = 1; i <= 6; i++) { for (j = 1; j <= 6; j++) { am[i, j] = 0.0; } } //要素長さ i = kakom[ne, 1]; j = kakom[ne, 2]; AL = Math.Sqrt((x[j] - x[i]) * (x[j] - x[i]) + (y[j] - y[i]) * (y[j] - y[i])); am[1, 1] = 1.0 / 3.0; am[2, 1] = 0.0; am[3, 1] = 0.0; am[4, 1] = 1.0 / 6.0; am[5, 1] = 0.0; am[6, 1] = 0.0; am[2, 2] = 13.0 / 35.0; am[3, 2] = 11.0 * AL / 210.0; am[4, 2] = 0.0; am[5, 2] = 9.0 / 70.0; am[6, 2] = -13.0 * AL / 420.0; am[3, 3] = AL * AL / 105.0; am[4, 3] = 0.0; am[5, 3] = 13.0 * AL / 420.0; am[6, 3] = -AL * AL / 140.0; am[4, 4] = 1.0 / 3.0; am[5, 4] = 0.0; am[6, 4] = 0.0; am[5, 5] = 13.0 / 35.0; am[6, 5] = -11.0 * AL / 210.0; am[6, 6] = AL * AL / 105.0; //対称項セット for (i = 1; i <= 5; i++) { for (j = i + 1; j <= 6; j++) { am[i, j] = am[j, i]; } } for (i = 1; i <= 6; i++) { for (j = 1; j <= 6; j++) { am[i, j] = gamma[ne] * AA[ne] * AL * am[i, j] / ggg; } } } //------------------------------------------------------------------------------ private void GLTMX(int ne, int[,] kakom, double[] x, double[] y, double[,] gt, double[,] hen) //座標変換行列と全体座標系への行列変換 { int i, j, k; double[,] w = new double[7, 7]; double AL, cs, sn; //行列初期化 for (i = 1; i <= 6; i++) { for (j = 1; j <= 6; j++) { hen[i, j] = 0.0; } } //要素長さ i = kakom[ne, 1]; j = kakom[ne, 2]; AL = Math.Sqrt((x[j] - x[i]) * (x[j] - x[i]) + (y[j] - y[i]) * (y[j] - y[i])); //座標変換行列 cs = (x[j] - x[i]) / AL; sn = (y[j] - y[i]) / AL; hen[1, 1] = cs; hen[1, 2] = sn; hen[2, 1] = -sn; hen[2, 2] = cs; hen[3, 3] = 1.0; for (i = 1; i <= 3; i++) { for (j = 1; j <= 3; j++) { hen[i + 3, j + 3] = hen[i, j]; } } //全体座標系への変換 //[w]=[T]//*[gt] for (i = 1; i <= 6; i++) { for (j = 1; j <= 6; j++) { w[i, j] = 0.0; for (k = 1; k <= 6; k++) { w[i, j] = w[i, j] + hen[k, i] * gt[k, j]; } } } //[w]*[T]=[T]//[gt][T]：[sm]は全体座標系での要素剛性行列 for (i = 1; i <= 6; i++) { for (j = 1; j <= 6; j++) { gt[i, j] = 0.0; for (k = 1; k <= 6; k++) { gt[i, j] = gt[i, j] + w[i, k] * hen[k, j]; } } } } //-------------------------------------------------------------------------------- private void MV_FRM(int ne, double[,] am, double ggg, double[] gkh, double[] gkv, double[] fevec) { //慣性力ベクトル int i, j; double[] work = new double[7]; work[1] = ggg * gkh[ne]; work[2] = ggg * gkv[ne]; work[3] = 0.0; work[4] = ggg * gkh[ne]; work[5] = ggg * gkv[ne]; work[6] = 0.0; //全体座標系での慣性力ベクトル for (i = 1; i <= 6; i++) { fevec[i] = 0.0; for (j = 1; j <= 6; j++) { fevec[i] = fevec[i] + am[i, j] * work[j]; } } } //-------------------------------------------------------------------------------- private void TV_FRM(int ne, int[,] kakom, double[] deltaT, double[] Em, double[] AA, double[] alpha, double[] fevec, double[,] hen) { //温度荷重相当ベクトル int i, j; double tempa; double[] work = new double[7]; i = kakom[ne, 1]; j = kakom[ne, 2]; tempa = 0.5 * (deltaT[i] + deltaT[j]); //要素座標系での温度荷重相当ベクトル work[1] = -1.0 * Em[ne] * AA[ne] * alpha[ne] * tempa; work[2] = 0.0; work[3] = 0.0; work[4] = 1.0 * Em[ne] * AA[ne] * alpha[ne] * tempa; work[5] = 0.0; work[6] = 0.0; //全体座標系での温度荷重相当ベクトル //（座標変換行列は転置=逆行列を要素座標系でのベクトルに乗じる） for (i = 1; i <= 6; i++) { fevec[i] = 0.0; for (j = 1; j <= 6; j++) { fevec[i] = fevec[i] + hen[j, i] * work[j]; } } } //-------------------------------------------------------------------------------- private void CALST_FRM(int ne, int[,] kakom, double[] dis, double[,] sm, double[,] hen, double[] reac, double[,] tsm, double[] Em, double[] AA, double[] alpha, double[] deltaT) { //断面力算定 int i, j, k; double[] wd = new double[7]; double[] ws = new double[7]; int ia, ir, iw; double tempa; //要素節点変位取得（全体座標系） for (i = 1; i <= 2; i++) { ia = kakom[ne, i]; wd[3 * i - 2] = dis[3 * ia - 2]; wd[3 * i - 1] = dis[3 * ia - 1]; wd[3 * i] = dis[3 * ia]; } //全体座標系での要素節点力算定 for (i = 1; i <= 6; i++) { ws[i] = 0.0; for (k = 1; k <= 6; k++) { ws[i] = ws[i] + sm[i, k] * wd[k]; } } //要素断面力算定：座標変換行列×全体座標系要素節点力 //今後使用しない全体剛性行列格納配列[tsm]に断面力を格納 for (i = 1; i <= 6; i++) { tsm[ne, i] = 0.0; for (k = 1; k <= 6; k++) { tsm[ne, i] = tsm[ne, i] + hen[i, k] * ws[k]; } } //軸力の温度応力補正 i = kakom[ne, 1]; j = kakom[ne, 2]; tempa = 0.5 * (deltaT[i] + deltaT[j]); tsm[ne, 1] = tsm[ne, 1] + Em[ne] * AA[ne] * alpha[ne] * tempa; tsm[ne, 4] = tsm[ne, 4] - Em[ne] * AA[ne] * alpha[ne] * tempa; //全体座標系での内力ベクトル //（座標変換行列は転置=逆行列を要素座標系でのベクトルに乗じる） for (i = 1; i <= 6; i++) { ws[i] = 0.0; for (j = 1; j <= 6; j++) { ws[i] = ws[i] + hen[j, i] * tsm[ne,j]; } } //要素節点力の全体系へのセット for (i = 1; i <= 2; i++) { ia = kakom[ne, i]; for (j = 1; j <= 3; j++) { ir = 3 * (ia - 1) + j; iw = 3 * (i - 1) + j; reac[ir] = reac[ir] + ws[iw]; } } } //------------------------------------------------------------------------------ private void ASMAT(int ne, int nod, int nhen, int[,] kakom, double[,] sm, double[,] tsm) { //剛性行列組立[nod：1要素節点数，nhen：1節点の自由度] int i, j, k, l, ki, kj, ks, js, kik, kjl, ksk, jsl; for (i = 1; i <= nod; i++) { for (j = 1; j <= nod; j++) { ki = (kakom[ne, i] - 1) * nhen; kj = (kakom[ne, j] - 1) * nhen; ks = (i - 1) * nhen; js = (j - 1) * nhen; for (k = 1; k <= nhen; k++) { for (l = 1; l <= nhen; l++) { kik = ki + k; kjl = kj + l; ksk = ks + k; jsl = js + l; tsm[kik, kjl] = tsm[kik, kjl] + sm[ksk, jsl]; } } } } } //------------------------------------------------------------------------------ private void ASVEC(int ne, int nod, int nhen, int[,] kakom, double[] fevec, double[] ftvec) { //荷重ベクトル組立[nod：1要素節点数，nhen：1節点の自由度] int i, k, ki, ks, kik, ksk; for (i = 1; i <= nod; i++) { ki = (kakom[ne, i] - 1) * nhen; ks = (i - 1) * nhen; for (k = 1; k <= nhen; k++) { kik = ki + k; ksk = ks + k; ftvec[kik] = ftvec[kik] + fevec[ksk]; } } } //------------------------------------------------------------------------------ private void CHBAND10(int n1, int m1, double[,] array) { //三角行列 int i, j, k, mm, mn; double z; array[1, 1] = Math.Sqrt(array[1, 1]); for (j = 2; j <= m1; j++) { array[1, j] = array[1, j] / array[1, 1]; } for (i = 2; i <= n1; i++) { if (n1 - i + 1 > m1) { mm = m1; } else { mm = n1 - i + 1; } for (j = 1; j <= mm; j++) { z = array[i, j]; if (j != m1) { if (m1 - j + 1 < i) { mn = m1 - j + 1; } else { mn = i; } for (k = 2; k <= mn; k++) { z = z - array[i - k + 1, k] * array[i - k + 1, j + k - 1]; } if (j > 1) { array[i, j] = z / array[i, 1]; } else { array[i, 1] = Math.Sqrt(z); } } else { array[i, j] = z / array[i, 1]; } } } } //------------------------------------------------------------------------------ private void CHBAND11(int n1, int m1, double[,] array, double[] vector) { //前進消去・後退代入 int i, j, mm; double z; vector[1] = vector[1] / array[1, 1]; for (i = 2; i <= n1; i++) { z = vector[i]; if (i > m1) { mm = m1; } else { mm = i; } for (j = 2; j <= mm; j++) { z = z - array[i - j + 1, j] * vector[i - j + 1]; } vector[i] = z / array[i, 1]; } vector[n1] = vector[n1] / array[n1, 1]; for (i = n1 - 1; i >= 1; i--) { z = vector[i]; for (j = 2; j <= m1; j++) { if (i + j - 1 > n1) break; z = z - array[i, j] * vector[i + j - 1]; } vector[i] = z / array[i, 1]; } } //------------------------------------------------------------------------------ private int IBAND(int mm, double[,] array, double fact0) { int ib, i, j; //----------------------------------------------------- //バンド幅計算 //----------------------------------------------------- ib = 1; for (i = 1; i <= mm - 1; i++) { for (j = mm; j >= i + 1; j--) { if (fact0 <= Math.Abs(array[i, j])) break; } if (ib <= j - i + 1) ib = j - i + 1; } if (mm < ib) ib = mm; return ib; } //------------------------------------------------------------------------------ private void BAND(int mm, int ib, double[,] array) { int i, j, k, mn; double[] work = new double[mm + 1]; //----------------------------------------------------- //フルマトリックスのバンド化記憶 //----------------------------------------------------- for (i = 2; i <= mm; i++) { if (i + ib - 1 < mm) { mn = ib; } else { mn = mm - i + 1; } for (j = 1; j <= ib; j++) { work[j] = 0.0; } for (k = 1; k <= mn; k++) { work[k] = array[i, i + k - 1]; } for (j = 1; j <= ib; j++) { array[i, j] = work[j]; } } } //------------------------------------------------------------------------------ private void IERROR2(string fnameW) { //結果出力（Appendモード） System.IO.StreamWriter sw; string dat; sw = new System.IO.StreamWriter(fnameW, true, System.Text.Encoding.Default); dat = "diagonal value is less or equal to 0"; sw.WriteLine(dat); sw.Close(); MessageBox.Show("対角項<=0", "異常終了通知"); this.Close(); } //------------------------------------------------------------------------------ } }