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 vcsSEEP3nod { 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() { //2次元飽和定常浸透流解析 int i, j, k, n, ia, ja, ne; System.IO.StreamReader sr; string dat; string[] sbuf; char delim = ','; string fnameR = ""; string fnameW = ""; string strcom; //書出用コメント int idan; //0:水平断面,1:鉛直断面 int NODT; //節点総数 int NELT; //要素総数 int MATEL; //材料種類数 int KOH; //全水頭指定節点数 int KOQ; //流量指定節点数 int nt; //全自由度(総節点数×1[1節点自由度]) int mm; //温度固定点処理後自由度 int ib; //バンド幅 double elarea; //要素面積 double fact0 = 1.0E-30; //0判定 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]; //------------------------------------------------------------------------------------------ //解析断面(0:水平,1:鉛直),節点数,要素数,材料種類数,全水頭固定節点数,流量指定節点数 //------------------------------------------------------------------------------------------ dat = sr.ReadLine(); sbuf = dat.Split(delim); idan = int.Parse(sbuf[0]); NODT = int.Parse(sbuf[1]); NELT = int.Parse(sbuf[2]); MATEL = int.Parse(sbuf[3]); KOH = int.Parse(sbuf[4]); KOQ = int.Parse(sbuf[5]); //------------------------------------------------------------------------------------------ //配列寸法宣言 //------------------------------------------------------------------------------------------ int[,] kakom = new int[NELT + 1, 4]; //要素構成接点番号 double[] Akx = new double[NELT + 1]; //要素x方向透水係数 double[] Aky = new double[NELT + 1]; //要素y方向透水係数 int[] matno = new int[NELT + 1]; //材料種別No double[,] wm = new double[MATEL + 1, 3]; //作業用材料物性記憶 double[] x = new double[NODT + 1]; //節点x座標 double[] y = new double[NODT + 1]; //節点y座標 int[] nokh = new int[KOH + 1]; //全水頭指定節点番号 int[] nokq = new int[KOQ + 1]; //流量指定節点番号 double[] Hinp = new double[KOH + 1]; //指定全水頭入力 double[] Qinp = new double[KOQ + 1]; //指定流量入力 double[] vx = new double[NELT + 1]; //x方向流速ベクトル double[] vy = new double[NELT + 1]; //y方向流速ベクトル nt = NODT * 1; int[] index = new int[nt + 1]; //温度固定点処理用 double[,] tak = new double[nt + 1, nt + 1]; //全体透水行列 double[] hvec = new double[nt + 1]; //全体全水頭ベクトル double[] qvec = new double[nt + 1]; //全体節点流量ベクトル double[] reac = new double[nt + 1]; //作業用ベクトル double[,] wak = new double[nt + 1, KOH + 1];//水頭指定節点対応列記憶 double[,] eak = new double[4, 4]; //要素透水行列 //------------------------------------------------------------------------------------------ //材料特性(Akx,Aky) //------------------------------------------------------------------------------------------ for (i = 1; i <= MATEL; i++) { dat = sr.ReadLine(); sbuf = dat.Split(delim); for (j = 1; j <= 2; j++) { wm[i, j] = double.Parse(sbuf[j - 1]); } } //------------------------------------------------------------------------------------------ //要素構成節点と材料種別(node-1,node-2,node-3,node-4,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]); kakom[ne, 3] = int.Parse(sbuf[2]); matno[ne] = int.Parse(sbuf[3]); for (i = 1; i <= MATEL; i++) { if (i == matno[ne]) { Akx[ne] = wm[i, 1]; Aky[ne] = wm[i, 2]; } } } //------------------------------------------------------------------------------------------ //節点座標(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 < KOH) { for (i = 1; i <= KOH; i++) { dat = sr.ReadLine(); sbuf = dat.Split(delim); nokh[i] = int.Parse(sbuf[0]); Hinp[i] = double.Parse(sbuf[1]); } } //------------------------------------------------------------------------------------------ //流量指定節点番号,指定流量 //------------------------------------------------------------------------------------------ if (0 < KOQ) { for (i = 1; i <= KOQ; i++) { dat = sr.ReadLine(); sbuf = dat.Split(delim); nokq[i] = int.Parse(sbuf[0]); Qinp[i] = double.Parse(sbuf[1]); } } sr.Close(); //----------------------------------------------------- //入力確認出力 //----------------------------------------------------- INPDAT(fnameW, strcom, idan, NODT, NELT, MATEL, KOH, KOQ, x, y, nokh, nokq, kakom, Akx, Aky, matno, Hinp, Qinp); //要素面積確認と異常終了通知 for (ne = 1; ne <= NELT; ne++) { i = kakom[ne, 1]; j = kakom[ne, 2]; k = kakom[ne, 3]; elarea = 0.5 * ((x[k] - x[j]) * y[i] + (x[i] - x[k]) * y[j] + (x[j] - x[i]) * y[k]); if (elarea < fact0) IERROR1(fnameW); } //***************************************************** //メイン処理計算 //***************************************************** //----------------------------------------------------- //全体行列クリア //----------------------------------------------------- for (i = 1; i <= nt; i++) { reac[i] = 0.0; hvec[i] = 0.0; qvec[i] = 0.0; for (j = 1; j <= nt; j++) { tak[i, j] = 0.0; } } //----------------------------------------------------- //全体行列組立 //----------------------------------------------------- for (ne = 1; ne <= NELT; ne++) { STIFF(ne, kakom, x, y, Akx, Aky, eak); ASMAT(ne, 3, 1, kakom, eak, tak); } //----------------------------------------------------- //全水頭指定節点対応列記憶,前処理 //----------------------------------------------------- if (0 < KOH) { for (k = 1; k <= KOH; k++) { for (i = 1; i <= nt; i++) { wak[i, k] = tak[i, nokh[k]]; } } for (k = 1; k <= KOH; k++) { for (i = 1; i <= nt; i++) { hvec[i] = hvec[i] + Hinp[k] * wak[i, k]; } } } if (0 < KOQ) { for (k = 1; k <= KOQ; k++) { qvec[nokq[k]] = Qinp[k]; } } for (i = 1; i <= nt; i++) { qvec[i] = qvec[i] - hvec[i]; } //----------------------------------------------------- //全水頭指定節点処理 //----------------------------------------------------- for (i = 1; i <= nt; i++) { index[i] = i; } if (0 < KOH) { for (i = 1; i <= KOH; i++) { n = nokh[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] = qvec[ia]; for (j = 1; j <= mm; j++) { ja = index[j]; tak[i, j] = tak[ia, ja]; } } //----------------------------------------------------- //バンド幅計算とフルマトリックスのバンド化記憶 //----------------------------------------------------- ib = IBAND(mm, tak, fact0); BAND(mm, ib, tak); //----------------------------------------------------- //対角項確認と異常終了通知 //----------------------------------------------------- for (i = 1; i <= mm; i++) { if (Math.Abs(tak[i, 1]) < fact0) IERROR2(fnameW); } //----------------------------------------------------- //コレスキー法 //----------------------------------------------------- CHBAND10(mm, ib, tak); //三角行列 CHBAND11(mm, ib, tak, reac); //----------------------------------------------------- //全水頭算定 //----------------------------------------------------- for (i = 1; i <= nt; i++) { hvec[i] = 0.0; } for (i = 1; i <= mm; i++) { ia = index[i]; hvec[ia] = reac[i]; } for (k = 1; k <= KOH; k++) { hvec[nokh[k]] = Hinp[k]; } //----------------------------------------------------- //節点流量算定 //----------------------------------------------------- for (i = 1; i <= nt; i++) { qvec[i] = 0.0; } if (0 < KOQ) { for (k = 1; k <= KOQ; k++) { qvec[nokq[k]] = Qinp[k]; } } if (0 < KOH) { for (k = 1; k <= KOH; k++) { for (i = 1; i <= nt; i++) { qvec[nokh[k]] = qvec[nokh[k]] + wak[i, k] * hvec[i]; } } } //----------------------------------------------------- //要素平均流速算定 //----------------------------------------------------- for (ne = 1; ne <= NELT; ne++) { CALV(ne, kakom, x, y, hvec, Akx, Aky, vx, vy); } //----------------------------------------------------- //時間計測完了 //----------------------------------------------------- entime = DateTime.Now; tspan = entime - sttime; label8.Text = entime.ToString(); label10.Text = "Calculation Time=" + tspan.TotalSeconds.ToString("0.00sec"); //***************************************************** //結果出力 //***************************************************** OUTDAT(fnameW, idan, NODT, NELT, x, y, hvec, qvec, vx, vy); label11.Text = "nt=" + nt.ToString() + " " + "mm=" + mm.ToString() + " " + "ib=" + ib.ToString(); MessageBox.Show("計算完了", "通知"); } //-------------------------------------------------------------------- private void INPDAT(string fnameW, string strcom, int idan, int NODT, int NELT, int MATEL, int KOH, int KOQ, double[] x, double[] y, int[] nokh, int[] nokq, int[,] kakom, double[] Akx, double[] Aky, int[] matno, double[] Hinp, double[] Qinp) { //入力確認出力 System.IO.StreamWriter sw; string dat; int i, ne; string[] w1 = new string[NODT + 1]; string[] w2 = new string[NODT + 1]; for (i = 1; i <= NODT; i++) { w1[i] = "--"; w2[i] = "--"; } if (0 < KOH) { for (i = 1; i <= KOH; i++) { w1[nokh[i]] = Hinp[i].ToString("E"); } } if (0 < KOQ) { for (i = 1; i <= KOQ; i++) { w2[nokq[i]] = Qinp[i].ToString("E"); } } sw = new System.IO.StreamWriter(fnameW, false, System.Text.Encoding.Default); dat = strcom; sw.WriteLine(dat); dat = "idan,NODT,NELT,MATEL,KOH,KOQ"; sw.WriteLine(dat); dat = idan.ToString("0"); dat = dat + "," + NODT.ToString("0"); dat = dat + "," + NELT.ToString("0"); dat = dat + "," + MATEL.ToString("0"); dat = dat + "," + KOH.ToString("0"); dat = dat + "," + KOQ.ToString("0"); sw.WriteLine(dat); dat = "*input-data"; sw.WriteLine(dat); dat = "node,x,y,Hinp,Qinp"; 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 + "," + w1[i]; dat = dat + "," + w2[i]; sw.WriteLine(dat); } dat = "element,node-1,node-2,node-3,Akx,Aky,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 + "," + kakom[ne, 3].ToString("0"); dat = dat + "," + Akx[ne].ToString("E"); dat = dat + "," + Aky[ne].ToString("E"); dat = dat + "," + matno[ne].ToString("0"); sw.WriteLine(dat); } sw.Close(); } //-------------------------------------------------------------------- private void OUTDAT(string fnameW, int idan, int NODT, int NELT, double[] x, double[] y, double[] hvec, double[] qvec, double[] vx, double[] vy) { //入力確認出力 System.IO.StreamWriter sw; string dat; int i, ne; sw = new System.IO.StreamWriter(fnameW, true, System.Text.Encoding.Default); dat = "*output-data"; sw.WriteLine(dat); dat = "node,x,y,hvec,qvec,P-head"; sw.WriteLine(dat); //節点番号,x座標,y座標,全水頭,節点流量,圧力水頭 //接点流量は,流入が正・流出が負の値 for (i = 1; i <= NODT; i++) { dat = i.ToString("0"); dat = dat + "," + x[i].ToString("E"); dat = dat + "," + y[i].ToString("E"); dat = dat + "," + hvec[i].ToString("E"); dat = dat + "," + qvec[i].ToString("E"); switch (idan) { case 0: dat = dat + "," + "--"; break; case 1: dat = dat + "," + (hvec[i] - y[i]).ToString("E"); break; } sw.WriteLine(dat); } dat = "element,vx,vy"; sw.WriteLine(dat); for (ne = 1; ne <= NELT; ne++) { dat = ne.ToString("0"); dat = dat + "," + vx[ne].ToString("E"); dat = dat + "," + vy[ne].ToString("E"); sw.WriteLine(dat); } dat = "Calculation time="+label10.Text; sw.WriteLine(dat); sw.Close(); } //-------------------------------------------------------------------- private void STIFF(int ne, int[,] kakom, double[] x, double[] y, double[] Akx, double[] Aky, double[,] eak) { //要素透水行列 int i,j,k; double b1,b2,b3,c1,c2,c3,Delta; double[,] wx=new double[4,4]; double[,] wy=new double[4,4]; for( i = 1;i<= 3;i++){ for(j = 1;j<=3;j++){ eak[i, j] = 0.0; wx[i, j] = 0.0; wy[i, j] = 0.0; } } i = kakom[ne, 1]; j = kakom[ne, 2]; k = kakom[ne, 3]; Delta = 0.5 * ((x[k] - x[j]) * y[i] + (x[i] - x[k]) * y[j] + (x[j] - x[i]) * y[k]); b1 = y[j] - y[k] ; c1 = x[k] - x[j]; b2 = y[k] - y[i] ; c2 = x[i] - x[k]; b3 = y[i] - y[j] ; c3 = x[j] - x[i]; wx[1, 1] = b1 * b1; wx[1, 2] = b1 * b2; wx[1, 3] = b1 * b3; wx[2, 1] = b2 * b1; wx[2, 2] = b2 * b2; wx[2, 3] = b2 * b3; wx[3, 1] = b3 * b1; wx[3, 2] = b3 * b2; wx[3, 3] = b3 * b3; wy[1, 1] = c1 * c1; wy[1, 2] = c1 * c2; wy[1, 3] = c1 * c3; wy[2, 1] = c2 * c1; wy[2, 2] = c2 * c2; wy[2, 3] = c2 * c3; wy[3, 1] = c3 * c1; wy[3, 2] = c3 * c2; wy[3, 3] = c3 * c3; for(i = 1;i<=3;i++){ for(j = 1;j<=3;j++){ eak[i, j] = Akx[ne] / 4.0 / Delta * wx[i, j] + Aky[ne] / 4.0 / Delta * wy[i, j]; } } } //-------------------------------------------------------------------- private void CALV(int ne, int[,] kakom, double[] x, double[] y, double[] hvec, double[] Akx, double[] Aky, double[] vx, double[] vy) { //流速ベクトル算定 int i,j,k; double b1,b2,b3,c1,c2,c3,Delta; double h1,h2,h3; i = kakom[ne, 1]; j = kakom[ne, 2]; k = kakom[ne, 3]; Delta = 0.5 * ((x[k] - x[j]) * y[i] + (x[i] - x[k]) * y[j] + (x[j] - x[i]) * y[k]); b1 = y[j] - y[k] ; c1 = x[k] - x[j]; b2 = y[k] - y[i] ; c2 = x[i] - x[k]; b3 = y[i] - y[j] ; c3 = x[j] - x[i]; h1 = hvec[i]; h2 = hvec[j]; h3 = hvec[k]; vx[ne] = -Akx[ne] / 2.0 / Delta * (b1 * h1 + b2 * h2 + b3 * h3); vy[ne] = -Aky[ne] / 2.0 / Delta * (c1 * h1 + c2 * h2 + c3 * h3); } //-------------------------------------------------------------------- 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 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 IERROR1(string fnameW) { //結果出力(Appendモード) System.IO.StreamWriter sw; string dat; sw = new System.IO.StreamWriter(fnameW, true, System.Text.Encoding.Default); dat = "area of element is less or equal to 0"; sw.WriteLine(dat); sw.Close(); MessageBox.Show("要素面積<=0", "異常終了通知"); this.Close(); } //------------------------------------------------------------------------------ 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(); } //------------------------------------------------------------------------------ } }