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 WindowsFormsApplication1 { public partial class Form1 : Form { public Form1() { InitializeComponent(); } private void toolStripButton1_Click(object sender, EventArgs e) { System.IO.StreamReader sr; string fnameL = ""; string fnameR = ""; string fnameW = ""; string dat; string[] sbuf; char delim = ' '; //入出力ファイル指定 openFileDialog1.InitialDirectory = System.IO.Directory.GetCurrentDirectory(); if (openFileDialog1.ShowDialog() == DialogResult.OK) { fnameL = openFileDialog1.FileName; } //データ入力 sr = new System.IO.StreamReader(fnameL, System.Text.Encoding.Default); while (!sr.EndOfStream) { dat = sr.ReadLine(); sbuf = dat.Split(delim); fnameR = System.IO.Path.GetDirectoryName(fnameL) + "\\" + sbuf[0]; fnameW = System.IO.Path.GetDirectoryName(fnameL) + "\\" + sbuf[1]; main_part(fnameR, fnameW); } sr.Close(); MessageBox.Show("計算完了", "通知"); } //-------------------------------------------------------------------------------- //-------------------------------------------------------------------------------- void main_part(string fnameR, string fnameW) { //平面骨組固有値解析 System.IO.StreamReader sr; System.IO.StreamWriter sw; string dat = ""; string[] sbuf; char delim = ','; int i, j, ne, n, ia, ja; 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 = 1.0e-30; //0判定基準 int jdg; //特異行列判定 DateTime sttime; //計算開始時刻 DateTime entime; //計算終了時刻 TimeSpan tspan; //計算時間 //時間計測開始 sttime = DateTime.Now; //***************************************************** //データ入力 //***************************************************** 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]; //要素断面二次モーメント int[] matno = new int[NELT + 1]; //材料種別No double[,] wm = new double[MATEL + 1, 4]; //材料種別作業用 double[] x = new double[NODT + 1]; //初期節点x座標 double[] y = new double[NODT + 1]; //初期節点y座標 int[] nokx = new int[NODT + 1]; //x方向変位既知節点番号 int[] noky = new int[NODT + 1]; //y方向変位既知節点番号 int[] nokz = new int[NODT + 1]; //回転変位既知節点番号 int[] index = new int[nt + 1]; //変位既知節点処理のためのインデックス double[,] tsm = new double[nt + 1, nt + 1]; //線形全体剛性マトリックス double[,] tsg = new double[nt + 1, nt + 1]; //非線形全体剛性マトリックス double[,] vec = new double[nt + 1, nt + 1]; //固有ベクトル double[] df = new double[nt + 1]; //全体節点外力ベクトル増分（入力値） double[] ftvec = new double[nt + 1]; //全体節点力ベクトル double[] dis = new double[nt + 1]; //全体節点変位増分 double[] reac = new double[nt + 1]; //全体節点変位増分 double[] ev = new double[nt + 1]; //固有値 double[,] sme = new double[7, 7]; //要素剛性マトリックス double[,] smg = new double[7, 7]; //要素剛性マトリックス double[] fevec = new double[7]; //要素節点荷重ベクトル double[,] hen = new double[7, 7]; //座標変換行列 double[,] reforce = new double[NELT + 1, 7]; //断面力 //------------------------------------------------------------------------------------------ //材料特性（Em,AA,AI） //------------------------------------------------------------------------------------------ for (i = 1; i <= MATEL; i++) { dat = sr.ReadLine(); sbuf = dat.Split(delim); wm[i, 1] = double.Parse(sbuf[0]); wm[i, 2] = double.Parse(sbuf[1]); wm[i, 3] = double.Parse(sbuf[2]); } //------------------------------------------------------------------------------------------ //要素構成節点と材料種別（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]; } } } //-------------------------------------------------------------------------------------- //節点座標（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]); } //-------------------------------------------------------------------------------------- //変位既知節点番号，既知節点変位量 //-------------------------------------------------------------------------------------- if (0 < KOX) { for (i = 1; i <= KOX; i++) { dat = sr.ReadLine(); sbuf = dat.Split(delim); nokx[i] = int.Parse(sbuf[0]); } } if (0 < KOY) { for (i = 1; i <= KOY; i++) { dat = sr.ReadLine(); sbuf = dat.Split(delim); noky[i] = int.Parse(sbuf[0]); } } if (0 < KOZ) { for (i = 1; i <= KOZ; i++) { dat = sr.ReadLine(); sbuf = dat.Split(delim); nokz[i] = int.Parse(sbuf[0]); } } //-------------------------------------------------------------------------------------- //節点番号,x方向荷重,y方向荷重,モーメント荷重 //-------------------------------------------------------------------------------------- for (i = 1; i <= nt; i++) { df[i] = 0.0; } if (0 < NF) { for (i = 1; i <= NF; i++) { dat = sr.ReadLine(); sbuf = dat.Split(delim); n = int.Parse(sbuf[0]); df[3 * n - 2] = double.Parse(sbuf[1]); df[3 * n - 1] = double.Parse(sbuf[2]); df[3 * n] = double.Parse(sbuf[3]); } } sr.Close(); //***************************************************** //入力確認出力 //***************************************************** INPDAT(fnameW, strcom, NODT, NELT, MATEL, KOX, KOY, KOZ, NF, x, y, df, nokx, noky, nokz, kakom, Em, AA, AI, matno); //***************************************************** //メイン処理計算 //***************************************************** //----------------------------------------------------- //変位既知節点自由度記憶 //----------------------------------------------------- for (i = 1; i <= nt; i++) { index[i] = 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; } } //***************************************************** //初期応力計算 {dis}=[K]^(-1){df} //***************************************************** for (ne = 1; ne <= NELT; ne++) { for (i = 1; i <= 6; i++) { reforce[ne, i] = 0.0; } } for (i = 1; i <= nt; i++) { for (j = 1; j <= nt; j++) { tsm[i, j] = 0.0; } } for (ne = 1; ne <= NELT; ne++) { FRAMESM(ne, kakom, x, y, Em, AA, AI, reforce, sme, smg); ASMAT(ne, 2, 3, kakom, sme, tsm); } for (i = 1; i <= nt; i++) { reac[i] = 0.0; dis[i] = 0.0; ftvec[i] = 0.0; } for (i = 1; i <= mm; i++) { ia = index[i]; reac[i] = df[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); jdg = CHBAND10(mm, ib, tsm);//三角行列 if (jdg == 1) goto epart; CHBAND11(mm, ib, tsm, reac);//前進消去・後退代入 for (i = 1; i <= mm; i++) { ia = index[i]; dis[ia] = reac[i]; } //内力計算 for (ne = 1; ne <= NELT; ne++) { CALST(ne, kakom, x, y, dis, Em, AA, AI, fevec, reforce); ASVEC(ne, 2, 3, kakom, fevec, ftvec); } DATAOUT_ISA(fnameW, NODT, NELT, x, y, dis, df, ftvec, reforce); //***************************************************** //固有値解析 [KL]{x}=ev[KG]{x} //***************************************************** for (i = 1; i <= nt; i++) { for (j = 1; j <= nt; j++) { tsm[i, j] = 0.0; tsg[i, j] = 0.0; } } for (ne = 1; ne <= NELT; ne++) { FRAMESM(ne, kakom, x, y, Em, AA, AI, reforce, sme, smg); ASMAT(ne, 2, 3, kakom, sme, tsm); ASMAT(ne, 2, 3, kakom, smg, tsg); } for (i = 1; i <= nt; i++) { ev[i] = 0.0; for (j = 1; j <= nt; j++) { vec[i, j] = 0.0; } } for (i = 1; i <= mm; i++) { ia = index[i]; for (j = 1; j <= mm; j++) { ja = index[j]; tsm[i, j] = tsm[ia, ja]; tsg[i, j] = tsg[ia, ja]; } } GJACOBI(mm, tsm, tsg, ev, vec); for (i = 1; i <= nt; i++) { dis[i] = 0.0; for (j = 1; j <= nt; j++) { tsm[i, j] = 0.0; } } for (i = 1; i <= mm; i++) { for (j = 1; j <= mm; j++) { ja = index[j]; tsm[ja, i] = vec[j, i]; } } DATAOUT_EVA(fnameW, NODT, NELT, x, y, ev, tsm); epart: //計算完了 entime = DateTime.Now; tspan = entime - sttime; sw = new System.IO.StreamWriter(fnameW, true, System.Text.Encoding.Default); dat = "NODT=" + NODT.ToString("0"); dat = dat + " nt=" + nt.ToString("0"); dat = dat + " mm=" + mm.ToString("0"); dat = dat + " ib=" + ib.ToString("0"); sw.WriteLine(dat); dat = "Calculation Time=" + tspan.TotalSeconds.ToString("0.00sec"); sw.WriteLine(dat); dat = entime.ToString(); sw.WriteLine(dat); sw.Close(); } //-------------------------------------------------------------------------------- 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[] df, int[] nokx, int[] noky, int[] nokz, int[,] kakom, double[] Em, double[] AA, double[] AI, 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 = "*input-data"; sw.WriteLine(dat); dat = "node,x,y,fx,fy,fz,x-fix,y-fix,z-fix"; 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 + "," + df[3 * i - 2].ToString("E"); dat = dat + "," + df[3 * i - 1].ToString("E"); dat = dat + "," + df[3 * i].ToString("E"); dat = dat + "," + k.ToString("0"); dat = dat + "," + l.ToString("0"); dat = dat + "," + m.ToString("0"); sw.WriteLine(dat); } dat = "element,node-1,node-2,E,A,I,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 + "," + matno[ne].ToString("0"); sw.WriteLine(dat); } sw.Close(); } //-------------------------------------------------------------------------------- void DATAOUT_ISA(string fnameW, int NODT, int NELT, double[] x, double[] y, double[] dis, double[] df, double[] ftvec, double[,] reforce) { //結果出力（Appendモード） int i, ne; System.IO.StreamWriter sw; string dat; sw = new System.IO.StreamWriter(fnameW, true, System.Text.Encoding.Default); dat = "***** Initial stress analysis *****"; sw.WriteLine(dat); dat = "*Displacements,forces"; sw.WriteLine(dat); dat = "node,x-coord,y-coord,x-dir,y-dir,z-dir,df-x,df-y,df-z,ft-x,ft-y,ft-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 + "," + df[3 * i - 2].ToString("E"); dat = dat + "," + df[3 * i - 1].ToString("E"); dat = dat + "," + df[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 resultant"; 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 + "," + reforce[ne, i].ToString("E"); } sw.WriteLine(dat); } sw.Close(); } //-------------------------------------------------------------------------------- void DATAOUT_EVA(string fnameW, int NODT, int NELT, double[] x, double[] y, double[] ev, double[,] vec) { System.IO.StreamWriter sw; string dat; int i, j, k, nt; nt = NODT * 3; sw = new System.IO.StreamWriter(fnameW, true, System.Text.Encoding.Default); dat = "***** Eigen value analysis *****"; sw.WriteLine(dat); i = 0; k = 0; do { i = i + 1; if (1e-30 < Math.Abs(ev[i]) && ev[i] < 0.0) { k = k + 1; dat = "*Eigen value"; sw.WriteLine(dat); dat = ev[i].ToString("E"); sw.WriteLine(dat); dat = "*Node data"; sw.WriteLine(dat); dat = "i,x,y,vec-x,vec-y,vec-z"; sw.WriteLine(dat); for (j = 1; j <= NODT; j++) { dat = j.ToString("0"); dat = dat + "," + x[j].ToString("E"); dat = dat + "," + y[j].ToString("E"); dat = dat + "," + vec[3 * j - 2, i].ToString("E"); dat = dat + "," + vec[3 * j - 1, i].ToString("E"); dat = dat + "," + vec[3 * j, i].ToString("E"); sw.WriteLine(dat); } } if (3 <= k) break; } while (i < nt); sw.Close(); } //-------------------------------------------------------------------------------- void FRAMESM(int ne, int[,] kakom, double[] x, double[] y, double[] Em, double[] AA, double[] AI, double[,] reforce, double[,] sme, double[,] smg) { int i, j, k; double[,] hen = new double[7, 7]; double[,] smu = new double[7, 7]; double[,] smv = new double[7, 7]; double[,] w = new double[7, 7]; double AL, P, BA, BI; double xx1, yy1, xx2, yy2; for (i = 1; i <= 6; i++) { for (j = 1; j <= 6; j++) { smv[i, j] = 0.0; smu[i, j] = 0.0; smg[i, j] = 0.0; sme[i, j] = 0.0; } } //軸力計算 i = kakom[ne, 1]; j = kakom[ne, 2]; xx1 = x[i]; yy1 = y[i]; xx2 = x[j]; yy2 = y[j]; AL = Math.Sqrt((xx2 - xx1) * (xx2 - xx1) + (yy2 - yy1) * (yy2 - yy1)); P = 0.5 * (reforce[ne, 4] - reforce[ne, 1]); //幾何剛性行列 smv[1, 1] = 0.0; smv[2, 1] = 0.0; smv[2, 2] = 6.0 / 5 / AL; smv[3, 1] = 0.0; smv[3, 2] = 1.0 / 10.0; smv[3, 3] = 2.0 * AL / 15.0; smv[4, 1] = 0.0; smv[4, 2] = 0.0; smv[4, 3] = 0.0; smv[4, 4] = 0.0; smv[5, 1] = 0.0; smv[5, 2] = -6.0 / 5.0 / AL; smv[5, 3] = -1.0 / 10.0; smv[5, 4] = 0.0; smv[5, 5] = 6.0 / 5.0 / AL; smv[6, 1] = 0.0; smv[6, 2] = 1.0 / 10.0; smv[6, 3] = -AL / 30.0; smv[6, 4] = 0.0; smv[6, 5] = -1.0 / 10.0; smv[6, 6] = 2.0 * AL / 15.0; for (i = 1; i <= 6; i++) { for (j = i; j <= 6; j++) { smv[i, j] = smv[j, i]; } } BA = AA[ne]; BI = AI[ne]; smu[1, 1] = 1.0 / AL; smu[2, 1] = 0.0; smu[2, 2] = 12.0 * BI / BA / AL / AL / AL; smu[3, 1] = 0.0; smu[3, 2] = 6 * BI / BA / AL / AL; smu[3, 3] = 4 * BI / BA / AL; smu[4, 1] = -1.0 / AL; smu[4, 2] = 0.0; smu[4, 3] = 0.0; smu[4, 4] = 1.0 / AL; smu[5, 1] = 0.0; smu[5, 2] = -12.0 * BI / BA / AL / AL / AL; smu[5, 3] = -6.0 * BI / BA / AL / AL; smu[5, 4] = 0.0; smu[5, 5] = 12.0 * BI / BA / AL / AL / AL; smu[6, 1] = 0.0; smu[6, 2] = 6.0 * BI / BA / AL / AL; smu[6, 3] = 2.0 * BI / BA / AL; smu[6, 4] = 0.0; smu[6, 5] = 6.0 * BI / BA / AL / AL; smu[6, 6] = 4.0 * BI / BA / AL; for (i = 1; i <= 6; i++) { for (j = i; j <= 6; j++) { smu[i, j] = smu[j, i]; } } for (i = 1; i <= 6; i++) { for (j = 1; j <= 6; j++) { smg[i, j] = smg[i, j] + P * smv[i, j]+P*smu[i,j]; } } HENMX(ne, xx1, yy1, xx2, yy2, hen); //[w]=[T]//*[k] 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] * smg[k, j]; } } } //[w]*[T]=[T]//[k,T]：[sm]は全体座標系での要素剛性行列 for (i = 1; i <= 6; i++) { for (j = 1; j <= 6; j++) { smg[i, j] = 0.0; for (k = 1; k <= 6; k++) { smg[i, j] = smg[i, j] + w[i, k] * hen[k, j]; } } } //線形弾性剛性行列 ELMSM(ne, xx1, yy1, xx2, yy2, Em, AA, AI, sme); //[w]=[T]//*[k] 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] * sme[k, j]; } } } //[w]*[T]=[T]//[k,T]：[sm]は全体座標系での要素剛性行列 for (i = 1; i <= 6; i++) { for (j = 1; j <= 6; j++) { sme[i, j] = 0.0; for (k = 1; k <= 6; k++) { sme[i, j] = sme[i, j] + w[i, k] * hen[k, j]; } } } } //-------------------------------------------------------------------------------- void ELMSM(int ne, double xx1, double yy1, double xx2, double yy2, double[] Em, double[] AA, double[] AI, double[,] sme) { //要素剛性行列 int i, j; double AL, EA, EI; AL = Math.Sqrt((xx2 - xx1) * (xx2 - xx1) + (yy2 - yy1) * (yy2 - yy1)); EA = Em[ne] * AA[ne]; EI = Em[ne] * AI[ne]; for (i = 1; i <= 6; i++) { for (j = 1; j <= 6; j++) { sme[i, j] = 0.0; } } sme[1, 1] = EA / AL; sme[2, 2] = 12.0 * EI / (AL * AL * AL); sme[3, 2] = 6.0 * EI / (AL * AL); sme[3, 3] = 4.0 * EI / AL; sme[4, 1] = -EA / AL; sme[4, 4] = sme[1, 1]; sme[5, 2] = -sme[2, 2]; sme[5, 3] = -sme[3, 2]; sme[5, 5] = sme[2, 2]; sme[6, 2] = sme[3, 2]; sme[6, 3] = 2.0 * EI / AL; sme[6, 5] = -sme[3, 2]; sme[6, 6] = sme[3, 3]; //対称項セット for (i = 1; i <= 6; i++) { for (j = 1; j <= i; j++) { sme[j, i] = sme[i, j]; } } } //-------------------------------------------------------------------------------- void CALST(int ne, int[,] kakom, double[] x, double[] y, double[] dis, double[] Em, double[] AA, double[] AI, double[] fevec, double[,] reforce) { //断面力算定 int i, k, ii, jj, ia; double[,] hen = new double[7, 7]; double[,] sme = new double[7, 7]; double[] wd = new double[7]; double[] wde = new double[7]; double[] ws = new double[7]; double xx1, yy1, xx2, yy2; ii = kakom[ne, 1]; jj = kakom[ne, 2]; xx1 = x[ii]; yy1 = y[ii]; xx2 = x[jj]; yy2 = y[jj]; HENMX(ne, xx1, yy1, xx2, yy2, hen); //要素節点変位（全体座標系） 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++) { wde[i] = 0.0; for (k = 1; k <= 6; k++) { wde[i] = wde[i] + hen[i, k] * wd[k]; } } ELMSM(ne, xx1, yy1, xx2, yy2, Em, AA, AI, sme); for (i = 1; i <= 6; i++) { ws[i] = 0.0; for (k = 1; k <= 6; k++) { ws[i] = ws[i] + sme[i, k] * wde[k]; } } //断面力保存 for (i = 1; i <= 6; i++) { reforce[ne, i] = ws[i]; } //全体座標系内力ベクトル for (i = 1; i <= 6; i++) { fevec[i] = 0.0; for (k = 1; k <= 6; k++) { fevec[i] = fevec[i] + hen[k, i] * reforce[ne, k]; } } } //-------------------------------------------------------------------------------- void HENMX(int ne, double xx1, double yy1, double xx2, double yy2, double[,] hen) { //座標変換行列 int i, j; double AL, cs, sn; AL = Math.Sqrt((xx2 - xx1) * (xx2 - xx1) + (yy2 - yy1) * (yy2 - yy1)); cs = (xx2 - xx1) / AL; sn = (yy2 - yy1) / AL; for (i = 1; i <= 6; i++) { for (j = 1; j <= 6; j++) { hen[i, j] = 0.0; } } 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]; } } } //-------------------------------------------------------------------------------- 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]; } } } } } //-------------------------------------------------------------------------------- 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]; } } } //-------------------------------------------------------------------------------- int 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 { if (z < 0) return 1; array[i, 1] = Math.Sqrt(z); } } else { array[i, j] = z / array[i, 1]; } } } return 0; } //-------------------------------------------------------------------------------- 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]; } } //-------------------------------------------------------------------------------- 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; } //-------------------------------------------------------------------------------- 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]; } } } //-------------------------------------------------------------------------------- void GJACOBI(int nnf, double[,] a, double[,] b, double[] ev, double[,] vec) { //************************************ // 一般化JACOBI法による固有値計算 // [パソコン構造力学，宮澤健二著] // 固有値ev[k],固有ベクトルvec[i,k] //************************************ int i, j, k, kmax, ip, iq, im; double fmk, fmg, epsk, epsg; double fact = 1.0E-20; double amax, dam, de1, de2, de3, dd, ss, x, gg, ww; fmk = 0.0; fmg = 0.0; for (i = 1; i <= nnf; i++) { for (j = 1; j <= nnf; j++) { if (Math.Abs(a[i, j]) > fmk) fmk = Math.Abs(a[i, j]); if (Math.Abs(b[i, j]) > fmg) fmg = Math.Abs(b[i, j]); } } epsk = fmk * fact; epsg = fmg * fact; for (i = 1; i <= nnf; i++) { for (j = 1; j <= nnf; j++) { vec[i, j] = 0.0; } vec[i, i] = 1.0; } kmax = 5 * nnf * nnf; for (k = 1; k <= kmax; k++) { ip = 1; iq = 2; amax = a[1, 2] * a[1, 2] / (a[1, 1] * a[1, 1] + a[2, 2] * a[2, 2]); for (i = 1; i <= nnf - 1; i++) { for (j = i + 1; j <= nnf; j++) { dam = a[i, j] * a[i, j] / (a[i, i] * a[i, i] + a[j, j] * a[j, j]); if (amax < dam) { ip = i; iq = j; amax = dam; } } } for (i = 1; i <= nnf - 1; i++) { for (j = i + 1; j <= nnf; j++) { dam = b[i, j] * b[i, j] / (b[i, i] * b[i, i] + b[j, j] * b[j, j]); if (amax < dam) { ip = i; iq = j; amax = dam; } } } de1 = a[ip, ip] * b[ip, iq] - a[ip, iq] * b[ip, ip]; de2 = a[iq, iq] * b[ip, iq] - a[ip, iq] * b[iq, iq]; de3 = a[ip, ip] * b[iq, iq] - a[iq, iq] * b[ip, ip]; dd = de3 * de3 + 4.0 * de1 * de2; ss = Math.Sqrt(dd); if (de3 < 0) { x = 0.5 * (de3 - ss); } else { x = 0.5 * (de3 + ss); } if (x == 0.0) { dam = 0.0; gg = -a[ip, iq] / a[iq, iq]; } else { dam = de2 / x; gg = -de1 / x; } for (j = 1; j <= nnf; j++) { ww = a[ip, j]; a[ip, j] = a[ip, j] + gg * a[iq, j]; a[iq, j] = a[iq, j] + dam * ww; ww = b[ip, j]; b[ip, j] = b[ip, j] + gg * b[iq, j]; b[iq, j] = b[iq, j] + dam * ww; } for (i = 1; i <= nnf; i++) { ww = a[i, ip]; a[i, ip] = a[i, ip] + gg * a[i, iq]; a[i, iq] = a[i, iq] + dam * ww; ww = b[i, ip]; b[i, ip] = b[i, ip] + gg * b[i, iq]; b[i, iq] = b[i, iq] + dam * ww; ww = vec[i, ip]; vec[i, ip] = vec[i, ip] + gg * vec[i, iq]; vec[i, iq] = vec[i, iq] + dam * ww; } if ((amax < epsk) && (amax < epsg)) break; } for (i = 1; i <= nnf; i++) { ev[i] = a[i, i] / b[i, i]; } //並び替え // ev[k]：最大固有値 k=1，最小固有値 k=nf // vec[i,k]：k番目の固有値に対応する固有ベクトル im = 0; for (k = 1; k <= nnf - 1; k++) { amax = -1.0E+30; for (i = k; i <= nnf; i++) { if (ev[i] >= amax) { im = i; amax = ev[i]; } } ww = ev[k]; ev[k] = ev[im]; ev[im] = ww; for (i = 1; i <= nnf; i++) { ww = vec[i, k]; vec[i, k] = vec[i, im]; vec[i, im] = ww; } } } //-------------------------------------------------------------------- } }