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 = ""; //個々出力ファイル int kprog; //第一中間点設定方法(1:方法1，2:方法2) double coefAB; //スケーリングパラメータ設定用係数 double deltaS0; //弧長 double coefS; //弧長に乗じる係数 int maxload; //荷重増分数最大値 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]; kprog = int.Parse(sbuf[2]); coefAB = double.Parse(sbuf[3]); deltaS0 = double.Parse(sbuf[4]); coefS = double.Parse(sbuf[5]); maxload = int.Parse(sbuf[6]); main_part(fnameR, fnameW, kprog,coefAB,deltaS0, coefS, maxload); } sr.Close(); MessageBox.Show("計算完了", "通知"); } //-------------------------------------------------------------------------------- int main_part(string fnameR, string fnameW, int kprog, double coefAB,double deltaS0, double coefS, int maxload) { //微小変形有限変位平面骨組構造解析 System.IO.StreamReader sr; System.IO.StreamWriter sw; string dat = ""; string[] sbuf; char delim = ','; int i, j, k, l, m, n, ne, 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 jflagc; //計算完了状態 int jdginv; //掃出法特異行列 double fend0, fend1; //荷重変化監視 double dend0, dend1; //変位変化監視 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[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[] df = new double[nt + 1]; //全体節点外力ベクトル増分（入力値） double[] ff = new double[nt + 1]; //全体節点外力ベクトル（荷重） double[] ftvec = new double[nt + 1]; //全体節点力ベクトル double[] dis = new double[nt + 1]; //全体節点変位増分 double[] dist = new double[nt + 1]; //全体節点変位 double[] dsav = new double[nt + 1]; //全体節点変位記憶 double[] dis0 = new double[nt + 1]; //全体節点変位増分 double[] disr = new double[nt + 1]; //全体節点変位増分 double[] work = new double[nt + 1]; //作業用ベクトル double[] reac = new double[nt + 1]; //作業用ベクトル double[,] sm = new double[7, 7]; //要素剛性マトリックス double[] fevec = new double[7]; //要素節点荷重ベクトル double[,] hen = new double[7, 7]; //座標変換行列 double dtest, dtest0; //変位収束判定用変数 double ftest, ftest0; //不平衡力収束判定用変数 int maxita = 50; //繰返数最大値 int itaration; //収束計算回数インクリメント int iload; //荷重増分インクリメント int nocon; //基準回数での収束か否かのフラグ double[,] reforce = new double[NELT + 1, 7]; //断面力 double[,] refsave = new double[NELT + 1, 7]; //釣合点断面力記憶 double lam, dlam; //弧長増分法パラメータ double deltaS; //弧長増分法パラメータ(計算に用いる弧長) double Spara; //弧長増分法パラメータ(スケーリングパラメータ) double fjdg0 = 1.0, fjdg; //弧長変更用 int jdg; //特異行列判定 //------------------------------------------------------------------------------------------ //材料特性（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(); //***************************************************** //入力確認出力 //***************************************************** 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(); //***************************************************** //メイン処理計算 //***************************************************** //----------------------------------------------------- //変位既知節点自由度記憶 //----------------------------------------------------- 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; } } //----------------------------------------------------- //ベクトルクリア //----------------------------------------------------- for (ne = 1; ne <= NELT; ne++) { for (i = 1; i <= 6; i++) { reforce[ne, i] = 0.0; refsave[ne, i] = 0.0; } } for (i = 1; i <= nt; i++) { ff[i] = 0.0; dist[i] = 0.0; dsav[i] = 0.0; dis0[i] = 0.0; disr[i] = 0.0; reac[i] = 0.0; } //----------------------------------------------------- //方法1スケーリングパラメータ設定 //----------------------------------------------------- ftest = 0.0; for (i = 1; i <= nt; i++) { ftest = ftest + df[i] * df[i]; } Spara = Math.Sqrt(coefAB / ftest); ftest = 0.0; //***************************************************** //荷重増分過程 //***************************************************** jflagc = 0; //正常終了フラグ初期設定 iload = 0; deltaS = deltaS0; ib = 0; fend0 = 0.0; jdginv = 0; do { iload = iload + 1; jdg = 0; again: itaration = 0; nocon = 0; lam = 0.0; dlam = 0.0; dtest0 = 0.0; ftest0 = 0.0; for (ne = 1; ne <= NELT; ne++) { for (i = 1; i <= 6; i++) { reforce[ne,i] = refsave[ne,i]; } } for (i = 1; i <= nt; i++) { dist[i] = dsav[i]; } /*初期処理 {df}=[KT]{dis0}*/ 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, dist, Em, AA, AI, reforce, sm); ASMAT(ne, 2, 3, kakom, sm, tsm); } for (i = 1; i <= nt; i++) { reac[i] = 0.0; dis0[i] = 0.0; work[i] = 0.0; } for (i = 1; i <= mm; i++) { ia = index[i]; work[i] = df[ia]; reac[i] = work[i]; for (j = 1; j <= mm; j++) { ja = index[j]; tsm[i,j] = tsm[ia,ja]; } } switch (jdg) { case 0: /*Cholesky*/ ib = IBAND(mm, tsm, fact0); BAND(mm, ib, tsm); jdg = CHBAND10(mm, ib, tsm);/*三角行列*/ CHBAND11(mm, ib, tsm, reac);/*前進消去・後退代入*/ break; case 1: /*Gauss-Jordan*/ jdginv = MATINV(mm, tsm); if (jdginv == 1) goto epart; for (i = 1; i <= mm; i++) { reac[i] = 0.0; for (j = 1; j <= mm; j++) { reac[i] = reac[i] + tsm[i,j] * work[j]; } } break; } for (i = 1; i <= mm; i++) { ia = index[i]; dis0[ia] = reac[i]; } switch(kprog){ case 1: dlam = CAL_LAM01(nt, deltaS, Spara, dis0, df); break; case 2: dlam = CAL_LAM02(nt, deltaS, coefAB,out Spara, dis0, df); break; } for (i = 1; i <= nt; i++) { dis[i] = dlam * dis0[i]; } /*====================================*/ /*収束過程*/ /*====================================*/ do { itaration = itaration + 1; lam = lam + dlam; for (i = 1; i <= nt; i++) { dist[i] = dist[i] + dis[i]; ftvec[i] = 0.0; } /*内力計算*/ for (ne = 1; ne <= NELT; ne++) { CALST(ne, kakom, x, y, dis, dist, Em, AA, AI, fevec, reforce); ASVEC(ne, 2, 3, kakom, fevec, ftvec); } /*不平衡力計算*/ for (i = 1; i <= nt; i++) { reac[i] = ff[i] + lam * df[i] - ftvec[i]; ftvec[i] = reac[i]; } for (i = 1; i <= mm; i++) { ia = index[i]; reac[i] = reac[ia]; dis[i] = dis[ia]; } /*収束判定*/ dtest = 0.0; ftest = 0.0; for (i = 1; i <= mm; i++) { dtest = dtest + Math.Abs(dis[i]); ftest = ftest + Math.Abs(reac[i]); } dtest0 = dtest0 + dtest; ftest0 = ftest0 + ftest; /*=====================================================*/ if (dtest < 1e-10) break; if (dtest < dtest0 * 1.0e-4 && ftest < ftest0 * 1.0e-4) break; /*=====================================================*/ for (i = 1; i <= nt; i++) { dis0[i] = 0.0; disr[i] = 0.0; reac[i] = 0.0; for (j = 1; j <= nt; j++) { tsm[i,j] = 0.0; } } for (ne = 1; ne <= NELT; ne++) { FRAMESM(ne, kakom, x, y, dist, Em, AA, AI, reforce, sm); ASMAT(ne, 2, 3, kakom, sm, tsm); } for (i = 1; i <= mm; i++) { ia = index[i]; reac[i] = df[ia]; work[i] = reac[i]; for (j = 1; j <= mm; j++) { ja = index[j]; tsm[i,j] = tsm[ia,ja]; } } switch (jdg) { case 0: /*Cholesky*/ ib = IBAND(mm, tsm, fact0); BAND(mm, ib, tsm); jdg = CHBAND10(mm, ib, tsm); if (jdg == 1) goto again; CHBAND11(mm, ib, tsm, reac);/*前進消去・後退代入*/ for (i = 1; i <= mm; i++) { ia = index[i]; dis0[ia] = reac[i]; reac[i] = ftvec[ia]; } CHBAND11(mm, ib, tsm, reac);/*前進消去・後退代入*/ break; case 1: /*Gauss-Jordan*/ jdginv = MATINV(mm, tsm); if (jdginv == 1) goto epart; for (i = 1; i <= mm; i++) { reac[i] = 0.0; for (j = 1; j <= mm; j++) { reac[i] = reac[i] + tsm[i,j] * work[j]; } } for (i = 1; i <= mm; i++) { ia = index[i]; dis0[ia] = reac[i]; work[i] = ftvec[ia]; } for (i = 1; i <= mm; i++) { reac[i] = 0.0; for (j = 1; j <= mm; j++) { reac[i] = reac[i] + tsm[i,j] * work[j]; } } break; } for (i = 1; i <= mm; i++) { ia = index[i]; disr[ia] = reac[i]; } dlam = CAL_LAMNEW(nt, Spara, dis0, disr, df); for (i = 1; i <= nt; i++) { dis[i] = dlam * dis0[i] + disr[i]; } } while (itaration <= maxita); /*====================================*/ if (maxita < itaration) nocon = 1; /*釣合点確定*/ for (i = 1; i <= nt; i++) { ff[i] = ff[i] + lam * df[i]; dsav[i] = dist[i]; } /*荷重変位確定*/ for (ne = 1; ne <= NELT; ne++) { for (i = 1; i <= 6; i++) { refsave[ne,i] = reforce[ne,i]; } } /*要素内力確定*/ DATAOUT(fnameW, NODT, NELT, iload, itaration, nocon, deltaS, lam, dtest, ftest, x, y, dist, ff, ftvec, reforce); /*荷重・変位が動かなくなったら計算終了*/ fend1 = 0.0; dend1 = 0.0; for (i = 1; i <= nt; i++) { fend1 = fend1 + Math.Abs(ff[i]); dend1 = dend1 + Math.Abs(dist[i]); } if (2 <= iload) { if (((1.0 - 1e-6) * fend0 < fend1 && fend1 < (1.0 + 1e-6) * fend0) && ((1.0 - 1e-6) * fend0 < fend1 && fend1 < (1.0 + 1e-6) * fend0)) { jflagc = 1;/*荷重・変位変化停止*/ goto epart; } } fend0 = fend1; dend0 = dend1; /*弧長変更判定*/ /*初回の増分荷重に対する増分荷重の比率が1/100未満となったら弧長を coefS 倍にして次回増分荷重を算定*/ if (iload == 1) { fjdg0 = 0.0; for (i = 1; i <= nt; i++) { if (fjdg0 < Math.Abs(lam * df[i]))fjdg0 = Math.Abs(lam * df[i]); } } fjdg = 0.0; for (i = 1; i <= nt; i++) { if (fjdg < Math.Abs(lam * df[i]))fjdg = Math.Abs(lam * df[i]); } if (fjdg < fjdg0 * 0.01) deltaS = coefS * deltaS; if (deltaS0 * 20.0 < deltaS) deltaS = deltaS0 * 20.0; } while (nocon == 0 && iload < maxload); /*******************************************************/ epart: if (nocon == 1) jflagc = 2;/*未収束*/ if (jdginv == 1) jflagc = 3;/*剛性行列特異*/ //計算完了 entime = DateTime.Now; tspan = entime - sttime; sw = new System.IO.StreamWriter(fnameW, true, System.Text.Encoding.Default); switch (jflagc) { case 0: dat = "jflagc=" + jflagc.ToString("0") + " : Normal completion"; break; case 1: dat = "jflagc=" + jflagc.ToString("0") + " : Stand still"; break; case 2: dat = "jflagc=" + jflagc.ToString("0") + " : No convergence"; break; case 3: dat = "jflagc=" + jflagc.ToString("0") + " : Singular matrix"; break; } sw.WriteLine(dat); 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 = "deltaS=" + deltaS.ToString("g"); dat = dat + " jdg=" + jdg.ToString(); dat = dat + " iload=" + iload.ToString(); dat = dat + " itaration=" + itaration.ToString(); sw.WriteLine(dat); dat = "Calculation Time=" + tspan.TotalSeconds.ToString("0.00sec"); sw.WriteLine(dat); dat = entime.ToString(); sw.WriteLine(dat); sw.Close(); return 0; } //-------------------------------------------------------------------------------- double CAL_LAM01(int nt, double deltaS, double Spara, double[] dis0, double[] df) { int i; double sum1, sum2; sum1 = 0.0; sum2 = 0.0; for (i = 1; i <= nt; i++) { sum1 = sum1 + dis0[i] * dis0[i]; sum2 = sum2 + df[i] * df[i]; } return Math.Sqrt(deltaS * deltaS / (sum1 + Spara * Spara * sum2)); } //-------------------------------------------------------------------------------- double CAL_LAM02(int nt, double deltaS, double coefAB,out double Spara, double[] dis0, double[] df) { int i; double sum1, sum2, sum3; sum1 = 0.0; sum2 = 0.0; sum3 = 0.0; for (i = 1; i <= nt; i++) { sum1 = sum1 + dis0[i] * dis0[i]; sum2 = sum2 + df[i] * df[i]; sum3 = sum3 + dis0[i] * df[i]; } Spara = 1.0 / coefAB * sum3 / sum2; return coefAB * Math.Sqrt(deltaS * deltaS / (sum1 + Spara * Spara * sum2)); } //-------------------------------------------------------------------------------- double CAL_LAMNEW(int nt, double Spara, double[] dis0, double[] disr, double[] df) { int i; double sum1, sum2, sum3; sum1 = 0.0; sum2 = 0.0; sum3 = 0.0; for (i = 1; i <= nt; i++) { sum1 = sum1 + dis0[i] * disr[i]; sum2 = sum2 + dis0[i] * dis0[i]; sum3 = sum3 + df[i] * df[i]; } return -sum1 / (sum2 + Spara * Spara * sum3); } //-------------------------------------------------------------------------------- void DATAOUT(string fnameW, int NODT, int NELT, int iload, int itaration, int nocon, double deltaS, double lam, double dtest, double ftest, double[] x, double[] y, double[] dist, double[] ff, 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); switch (nocon) { case 0: dat = "***** Output *****"; break; case 1: dat = "##### Output #####"; break; } sw.WriteLine(dat); dat = "iload,itaration,nocon,deltaS,lam,dtest,ftest"; sw.WriteLine(dat); dat = iload.ToString("0"); dat = dat + "," + itaration.ToString("0"); dat = dat + "," + nocon.ToString("0"); dat = dat + "," + deltaS.ToString("E"); dat = dat + "," + lam.ToString("E"); dat = dat + "," + dtest.ToString("E"); dat = dat + "," + ftest.ToString("E"); sw.WriteLine(dat); dat = "*Node values"; sw.WriteLine(dat); dat = "node,x-coord,y-coord,x-dir,y-dir,z-dir,ff-x,ff-y,ff-z,ftvec-x,ftvec-y,ftvec-z"; sw.WriteLine(dat); for (i = 1; i <= NODT; i++) { dat = i.ToString("0"); dat = dat + "," + (x[i] + dist[3 * i - 2]).ToString("E"); dat = dat + "," + (y[i] + dist[3 * i - 1]).ToString("E"); dat = dat + "," + dist[3 * i - 2].ToString("E"); dat = dat + "," + dist[3 * i - 1].ToString("E"); dat = dat + "," + dist[3 * i].ToString("E"); dat = dat + "," + ff[3 * i - 2].ToString("E"); dat = dat + "," + ff[3 * i - 1].ToString("E"); dat = dat + "," + ff[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 FRAMESM(int ne, int[,] kakom, double[] x, double[] y, double[] dist, double[] Em, double[] AA, double[] AI, double[,] reforce, double[,] sm) { int i, j, k; double[,] hen = new double[7, 7]; double[,] sme = new double[7, 7]; double[,] smv = new double[7, 7]; double[,] w = new double[7, 7]; double AL, P; double xx1, yy1, xx2, yy2; for (i = 1; i <= 6; i++) { for (j = 1; j <= 6; j++) { smv[i, j] = 0.0; sm[i, j] = 0.0; } } //軸力計算 i = kakom[ne, 1]; j = kakom[ne, 2]; xx1 = x[i] + dist[3 * i - 2]; yy1 = y[i] + dist[3 * i - 1]; xx2 = x[j] + dist[3 * j - 2]; yy2 = y[j] + dist[3 * j - 1]; 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]; } } //線形弾性剛性行列 ELMSM(ne, xx1, yy1, xx2, yy2, Em, AA, AI, sme); for (i = 1; i <= 6; i++) { for (j = 1; j <= 6; j++) { sm[i, j] = sme[i, j] + P * smv[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] * sm[k, j]; } } } //[w]*[T]=[T]//[k,T]：[sm]は全体座標系での要素剛性行列 for (i = 1; i <= 6; i++) { for (j = 1; j <= 6; j++) { sm[i, j] = 0.0; for (k = 1; k <= 6; k++) { sm[i, j] = sm[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[] dist, 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 AL, BL, deltai, deltaj; double[] wd = new double[7]; double[] wde = new double[7]; double[] ws = new double[7]; double xx1, yy1, xx2, yy2; ROTATE(ne, kakom, x, y, dis, dist, out deltai, out deltaj); //要素節点変位増分（全体座標系） 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]; } ii = kakom[ne, 1]; jj = kakom[ne, 2]; //前回要素座標系からの変位増分 xx1 = x[ii] + dist[3 * ii - 2] - dis[3 * ii - 2]; yy1 = y[ii] + dist[3 * ii - 1] - dis[3 * ii - 1]; xx2 = x[jj] + dist[3 * jj - 2] - dis[3 * jj - 2]; yy2 = y[jj] + dist[3 * jj - 1] - dis[3 * jj - 1]; HENMX(ne, xx1, yy1, xx2, yy2, hen); 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]; } } //前回部材長 xx1 = x[ii] + dist[3 * ii - 2] - dis[3 * ii - 2]; yy1 = y[ii] + dist[3 * ii - 1] - dis[3 * ii - 1]; xx2 = x[jj] + dist[3 * jj - 2] - dis[3 * jj - 2]; yy2 = y[jj] + dist[3 * jj - 1] - dis[3 * jj - 1]; BL = Math.Sqrt((xx2 - xx1) * (xx2 - xx1) + (yy2 - yy1) * (yy2 - yy1)); //今回部材長 xx1 = x[ii] + dist[3 * ii - 2]; yy1 = y[ii] + dist[3 * ii - 1]; xx2 = x[jj] + dist[3 * jj - 2]; yy2 = y[jj] + dist[3 * jj - 1]; AL = Math.Sqrt((xx2 - xx1) * (xx2 - xx1) + (yy2 - yy1) * (yy2 - yy1)); wd[1] = 0.0; wd[2] = 0.0; wd[3] = deltai; wd[4] = AL - BL; wd[5] = 0.0; wd[6] = deltaj; //要素内力増分(前回微小変形要素剛性使用) xx1 = x[ii] + dist[3 * ii - 2] - dis[3 * ii - 2]; yy1 = y[ii] + dist[3 * ii - 1] - dis[3 * ii - 1]; xx2 = x[jj] + dist[3 * jj - 2] - dis[3 * jj - 2]; yy2 = y[jj] + dist[3 * jj - 1] - dis[3 * jj - 1]; 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] * wd[k]; } } //断面力保存 for (i = 1; i <= 6; i++) { reforce[ne, i] = reforce[ne, i] + ws[i]; } //全体座標系での変形後要素内力 xx1 = x[ii] + dist[3 * ii - 2]; yy1 = y[ii] + dist[3 * ii - 1]; xx2 = x[jj] + dist[3 * jj - 2]; yy2 = y[jj] + dist[3 * jj - 1]; HENMX(ne, xx1, yy1, xx2, yy2, hen); 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 ROTATE(int ne, int[,] kakom, double[] x, double[] y, double[] dis, double[] dist, out double deltai, out double deltaj) { //剛体回転除去 int i, k, ii, jj, ia; double[,] hen = new double[7, 7]; double[] wd = new double[7]; double[] wde = new double[7]; double xx1, yy1, xx2, yy2; double BL, Ti, Tj, theta1i, theta1j, theta2i, theta2j; ii = kakom[ne, 1]; jj = kakom[ne, 2]; //前回における初期部材座標系からの変位量 for (i = 1; i <= 2; i++) { ia = kakom[ne, i]; wd[3 * i - 2] = dist[3 * ia - 2] - dis[3 * ia - 2]; wd[3 * i - 1] = dist[3 * ia - 1] - dis[3 * ia - 1]; wd[3 * i] = dist[3 * ia] - dis[3 * ia]; } //初期部材長さ xx1 = x[ii]; yy1 = y[ii]; xx2 = x[jj]; yy2 = y[jj]; BL = Math.Sqrt((xx2 - xx1) * (xx2 - xx1) + (yy2 - yy1) * (yy2 - yy1)); HENMX(ne, xx1, yy1, xx2, yy2, hen); 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]; } } //前回の内力計算用回転角 Ti = Math.Tan(wde[3]); Tj = Math.Tan(wde[6]); theta1i = ((BL + wde[4] - wde[1]) * Ti - (wde[5] - wde[2])) / ((BL + wde[4] - wde[1]) + (wde[5] - wde[2]) * Ti); theta1j = ((BL + wde[4] - wde[1]) * Tj - (wde[5] - wde[2])) / ((BL + wde[4] - wde[1]) + (wde[5] - wde[2]) * Tj); //今回における初期部材座標系からの変位量 for (i = 1; i <= 2; i++) { ia = kakom[ne, i]; wd[3 * i - 2] = dist[3 * ia - 2]; wd[3 * i - 1] = dist[3 * ia - 1]; wd[3 * i] = dist[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]; } } //前回の内力計算用回転角 Ti = Math.Tan(wde[3]); Tj = Math.Tan(wde[6]); theta2i = ((BL + wde[4] - wde[1]) * Ti - (wde[5] - wde[2])) / ((BL + wde[4] - wde[1]) + (wde[5] - wde[2]) * Ti); theta2j = ((BL + wde[4] - wde[1]) * Tj - (wde[5] - wde[2])) / ((BL + wde[4] - wde[1]) + (wde[5] - wde[2]) * Tj); //内力計算用回転角増分 deltai = theta2i - theta1i; deltaj = theta2j - theta1j; } //-------------------------------------------------------------------------------- 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]; } } } //-------------------------------------------------------------------------------- int MATINV(int n, double[,] ary) { /*****************************************/ /* 逆行列計算 */ /* 山田嘉昭著：マトリックス法材料力学 */ //****************************************/ int i, j, k, l, l1, irow = 0, icol = 0, jrow = 0, jcol = 0; double amax, s, t; double determ = 1.0; double factor = 1.0e-30; int[] ipivot = new int[n + 1]; double[] pivot = new double[n + 1]; int[,] index = new int[n + 1, 3]; for (j = 1; j <= n; j++) { ipivot[j] = 0; } for (i = 1; i <= n; i++) { amax = 0.0; for (j = 1; j <= n; j++) { if (ipivot[j] - 1 != 0) { for (k = 1; k <= n; k++) { if (0 < ipivot[k] - 1) return 0; if (ipivot[k] - 1 < 0) { if (Math.Abs(amax) - Math.Abs(ary[j, k]) < 0) { irow = j; icol = k; amax = ary[j, k]; } } } } } ipivot[icol] = ipivot[icol] + 1; if ((irow + icol) % 2 != 0) { determ = -determ; } if (irow - icol != 0) { for (l = 1; l <= n; l++) { s = ary[irow, l]; ary[irow, l] = ary[icol, l]; ary[icol, l] = s; } } index[i, 1] = irow; index[i, 2] = icol; pivot[i] = ary[icol, icol]; if (Math.Abs(pivot[i]) < factor) return 1; if (determ <= 1.0E+30) { determ = determ * pivot[i]; } ary[icol, icol] = 1.0; for (l = 1; l <= n; l++) { ary[icol, l] = ary[icol, l] / pivot[i]; } for (l1 = 1; l1 <= n; l1++) { if (l1 - icol != 0) { t = ary[l1, icol]; ary[l1, icol] = 0.0; for (l = 1; l <= n; l++) { ary[l1, l] = ary[l1, l] - ary[icol, l] * t; } } } } for (i = 1; i <= n; i++) { l = n + 1 - i; if (index[l, 1] - index[l, 2] != 0) { jrow = index[l, 1]; jcol = index[l, 2]; for (k = 1; k <= n; k++) { s = ary[k, jrow]; ary[k, jrow] = ary[k, jcol]; ary[k, jcol] = s; } } } return 0; } /*--------------------------------------------------------------------------------*/ } }