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 vcsFEMondoCH
{
    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,k,n,nnn,ia,ja,ne,it;
            System.IO.StreamReader sr;
            System.IO.StreamWriter sw;
            string dat;
            string[] sbuf;
            char delim = ',';
            string fnameR="";
            string fnameW="";
            string ftinp="";

            string strcom;          //書出用コメント
            int NODT;               //節点総数
            int NELT;               //要素総数
            int MATEL;              //材料種類数
            int KOT;                //温度指定節点数
            int KOC;                //熱伝達指定辺数
            int nt;                 //全自由度（総節点数×１[１節点自由度]）
            int mm;                 //温度固定点処理後自由度
            int ib;                 //バンド幅
            double delta;           //時間刻み
            double ttime;           //時刻
            double elarea1;         //要素半面積
            double elarea2;         //要素半面積
            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];
            //--------------------------------------------------------------------------------------
            //時刻歴読み込みファイル指定
            //--------------------------------------------------------------------------------------
            dat = sr.ReadLine() ; sbuf = dat.Split(delim);
            ftinp = System.IO.Path.GetDirectoryName(fnameR) + "\\" + sbuf[0];
            //------------------------------------------------------------------------------------------
            //節点数，要素数，材料種類数,温度固定点数，熱伝達境界辺数，時間刻み
            //------------------------------------------------------------------------------------------
            dat = sr.ReadLine() ; sbuf = dat.Split(delim);
            NODT  = int.Parse(sbuf[0]);
            NELT  = int.Parse(sbuf[1]);
            MATEL = int.Parse(sbuf[2]);
            KOT   = int.Parse(sbuf[3]);
            KOC   = int.Parse(sbuf[4]);
            delta = double.Parse(sbuf[5]);
            //------------------------------------------------------------------------------------------
            //配列寸法宣言
            //------------------------------------------------------------------------------------------
            double alpc;   //作業用熱伝達率
            double tc;     //作業用熱伝達境界温度
            int kchen;     //作業用熱伝達境界辺番号
            double dotq;   //作業用発熱率
            int[,] kakom=new int[NELT+1, 5];     //要素構成節点番号
            double[] Ak=new double[NELT+1];      //要素熱伝導率
            double[] Ac=new double[NELT+1];      //要素比熱
            double[] Arho=new double[NELT+1];    //要素単位体積重量
            double[] Tk=new double[NELT+1];      //要素最終温度上昇量
            double[] Al=new double[NELT+1];      //要素発熱速度係数
            int[] matno=new int[NELT+1];         //材料種別No
            double[,] wm=new double[MATEL+1, 6]; //作業用材料物性記憶
            double[] x=new double[NODT+1];       //節点x座標
            double[] y=new double[NODT+1];       //節点y座標
            int[] nokt=new int[KOT+1];           //温度指定節点番号
            int[,] nekc=new int[KOC+1, 3];       //熱伝達境界要素・辺番号
            double[] alphac=new double[KOC+1];   //指定辺熱伝達率
            double[] Tinp=new double[KOT+1];     //指定温度入力
            double[] Tcin1=new double[KOC+1];    //指定熱伝達境界温度入力（前回）
            double[] Tcin2=new double[KOC+1];    //指定熱伝達境界温度入力（現在）
            nt = NODT * 1;
            int[] index=new int[nt+1];               //温度固定点処理用
            double[,] tck1=new double[nt+1, nt+1];   //全体剛性行列（前回）
            double[,] tck2=new double[nt+1, nt+1];   //全体剛性行列（現在）
            double[,] rtck2=new double[nt+1, KOT+1]; //全体剛性行列（現在）拘束節点列保存
            double[] tempe0 = new double[nt + 1];    //初期節点温度ベクトル
            double[] tempe = new double[nt + 1];     //全体節点温度ベクトル
            double[] reac=new double[nt+1];          //作業用節点温度ベクトル
            double[] ftvec=new double[nt+1];         //全体熱流束ベクトル

            int itmax;                      //時刻歴ステップ数
            int nout;                       //出力節点数
            int ntout;                      //全節点温度出力回数
            double[,] ck1=new double[5, 5]; //要素剛性行列（前回）
            double[,] ck2=new double[5, 5]; //要素剛性行列（現在）
            double[,] ht=new double[5, 5];  //要素熱伝達行列
            double[] fq=new double[5];      //要素発熱率ベクトル
            double[] fc=new double[5];      //要素熱伝達ベクトル
            //------------------------------------------------------------------------------------------
            //材料特性（Ak,Ac,Arho,Tk,Al）
            //------------------------------------------------------------------------------------------
            for(i=1;i<=MATEL;i++){
                dat = sr.ReadLine() ; sbuf = dat.Split(delim);
                for(j=1;j<=5;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]);
                kakom[ne, 4] = int.Parse(sbuf[3]);
                matno[ne]    = int.Parse(sbuf[4]);
                for(i=1;i<=MATEL;i++){
                    if(i==matno[ne]){
                        Ak[ne]   = wm[i, 1];
                        Ac[ne]   = wm[i, 2];
                        Arho[ne] = wm[i, 3];
                        Tk[ne]   = wm[i, 4];
                        Al[ne]   = wm[i, 5];
                    }
            }}
            //------------------------------------------------------------------------------------------
            //節点座標（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]);
                tempe0[i] = double.Parse(sbuf[2]);
            }
            //------------------------------------------------------------------------------------------
            //温度指定節点番号
            //------------------------------------------------------------------------------------------
            if(0<KOT){
                for(i=1;i<=KOT;i++){
                    dat = sr.ReadLine() ; sbuf = dat.Split(delim) ; nokt[i] = int.Parse(sbuf[0]);
                }
            }
            //------------------------------------------------------------------------------------------
            //熱伝達境界指定要素・辺番号・熱伝達率
            //------------------------------------------------------------------------------------------
            if(0<KOC){
                for(ne=1;ne<=KOC;ne++){
                    dat = sr.ReadLine() ; sbuf = dat.Split(delim);
                    nekc[ne, 1] = int.Parse(sbuf[0]); //要素番号
                    nekc[ne, 2] = int.Parse(sbuf[1]); //辺番号
                    alphac[ne]  = double.Parse(sbuf[2]); //熱伝達率
                }
            }
            //------------------------------------------------------------------------------------------
            //出力節点数,出力節点番号
            //------------------------------------------------------------------------------------------
            dat = sr.ReadLine() ; sbuf = dat.Split(delim);
            nout = int.Parse(sbuf[0]);
            int[] noout=new int[nout+1];//出力節点番号
            dat = sr.ReadLine() ; sbuf = dat.Split(delim);
            for(i=1;i<=nout;i++){
                noout[i] = int.Parse(sbuf[i - 1]);
            }
            //------------------------------------------------------------------------------------------
            //全節点温度出力回数,出力ステップ番号
            //------------------------------------------------------------------------------------------
            dat = sr.ReadLine() ; sbuf = dat.Split(delim);
            ntout = int.Parse(sbuf[0]);
            int[] nntout=new int[ntout+1];      //全節点温度出力ステップ
            double[,] tempend=new double[NODT+1,ntout+1];   //出力用節点温度
            if(0<ntout){
                dat = sr.ReadLine() ; sbuf = dat.Split(delim);
                for(i=1;i<=ntout;i++){
                    nntout[i] = int.Parse(sbuf[i - 1]);
                }
            }
            sr.Close();

            //---------------------------------------------------
            //時刻歴入力ファイル空読み
            //---------------------------------------------------
            sr = new System.IO.StreamReader(ftinp, System.Text.Encoding.Default);
            dat = sr.ReadLine(); //コメント
            itmax = 0;
            while (!sr.EndOfStream)
            {
                dat = sr.ReadLine(); sbuf = dat.Split(delim);
                itmax = int.Parse(sbuf[0]);
            }
            sr.Close();
            //*****************************************************
            //入力確認出力
            //*****************************************************
            INPDAT(fnameW, strcom, NODT, NELT, MATEL, KOT, KOC, delta,nout,ntout,itmax,
                   x, y, tempe, nokt, nekc, kakom, Ak, Ac, Arho, Tk, Al, matno);
            
            //要素面積確認と異常終了通知
            for (ne = 1; ne <= NELT; ne++)
            {
                i = kakom[ne, 1]; j = kakom[ne, 2]; k = kakom[ne, 3];
                elarea1 = 0.5 * ((x[k] - x[j]) * y[i] + (x[i] - x[k]) * y[j] + (x[j] - x[i]) * y[k]);
                if (elarea1 < fact0) IERROR1(fnameW);
                i = kakom[ne, 1]; j = kakom[ne, 3]; k = kakom[ne, 4];
                elarea2 = 0.5 * ((x[k] - x[j]) * y[i] + (x[i] - x[k]) * y[j] + (x[j] - x[i]) * y[k]);
                if (elarea2 < fact0) IERROR1(fnameW);
            }
            //*****************************************************
            //メイン処理計算
            //*****************************************************
            //-----------------------------------------------------
            //全体行列クリア
            //-----------------------------------------------------
            for(i=1;i<=nt;i++){
                for(j=1;j<=nt;j++){
                    tck1[i, j] = 0.0;
                    tck2[i, j] = 0.0;
            }}
            //-----------------------------------------------------
            //全体行列組立
            //-----------------------------------------------------
            for(ne=1;ne<=NELT;ne++){
                STIFF(ne, kakom, x, y, Ak, Ac, Arho, ck1, ck2, delta);
                ASMAT(ne, 4, 1, kakom, ck1, tck1); //前回剛性-0.5*[K]+1/dt*[C]
                ASMAT(ne, 4, 1, kakom, ck2, tck2); //現在剛性+0.5*[K]+1/dt*[C]
                //熱伝達境界処理
                if(0<KOC){
                    for(k=1;k<=KOC;k++){
                        if(ne==nekc[k, 1]){
                            kchen = nekc[k, 2];
                            alpc = alphac[k];
                            HTMAT(ne, kchen, kakom, x, y, alpc, ht); //熱伝達行列
                            for(i=1;i<=4;i++){
                                for(j=1;j<=4;j++){
                                    ht[i, j] = 0.5 * ht[i, j];
                            }}
                            ASMAT(ne, 4, 1, kakom, ht, tck2); //現在剛性
                            for(i=1;i<=4;i++){
                                for(j=1;j<=4;j++){
                                    ht[i, j] = -ht[i, j];
                            }}
                            ASMAT(ne, 4, 1, kakom, ht, tck1); //前回剛性
                        }
                    }
                }   
            }
            //-----------------------------------------------------
            //全体行列温度指定列要素保存
            //-----------------------------------------------------
            if(0<KOT){
                for(k=1;k<=KOT;k++){
                    for(i=1;i<=nt;i++){
                        rtck2[i, k] = tck2[i, nokt[k]];
                }}
            }
            //-----------------------------------------------------
            //温度指定節点処理
            //-----------------------------------------------------
            for(i=1;i<=nt;i++){index[i] = i ;}
            if(0<KOT){
                for(i=1;i<=KOT;i++){
                    n = nokt[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];
                for(j=1;j<=mm;j++){
                    ja = index[j];
                    tck2[i, j] = tck2[ia, ja];
            }}
            //-----------------------------------------------------
            //バンド幅計算とフルマトリックスのバンド化記憶
            //-----------------------------------------------------
            ib=IBAND(mm, tck2, fact0);
            BAND(mm, ib, tck2);
            //-----------------------------------------------------
            //対角項確認と異常終了通知
            //-----------------------------------------------------
            for (i = 1; i <= mm; i++)
            {
                if (Math.Abs(tck2[i, 1]) < fact0) IERROR2(fnameW);
            }
            //-----------------------------------------------------
            //コレスキー法（三角行列）
            //-----------------------------------------------------
            CHBAND10(mm, ib, tck2);       //三角行列
               
            //-----------------------------------------------------
            //時刻歴出力準備
            //-----------------------------------------------------
            /*初期温度ベクトルセット*/
            for (i = 1; i <= nt; i++)
            {
                tempe[i] = tempe0[i];
            }
            /*熱伝達境界初期値設定*/
            ia = 0; ja = 0;
            if (0 < KOC)
            {
                for (k = 1; k <= KOC; k++)
                {
                    ne = nekc[k,1];
                    switch (nekc[k,2])
                    {
                        case 1: ia = kakom[ne,1]; ja = kakom[ne,2]; break;
                        case 2: ia = kakom[ne,2]; ja = kakom[ne,3]; break;
                        case 3: ia = kakom[ne,3]; ja = kakom[ne,4]; break;
                        case 4: ia = kakom[ne,4]; ja = kakom[ne,1]; break;
                    }
                    Tcin1[k] = 0.5 * (tempe[ia] + tempe[ja]);
                }
            }
            
            sw = new System.IO.StreamWriter(fnameW, true, System.Text.Encoding.Default);
            dat = "*temperature records of selected nodes" ; sw.WriteLine(dat);
            //======================
            it = 0; ttime = 0.0;
            //======================
            dat = "it,time";
            for(i=1;i<=nout;i++){
                dat = dat + ",Node-" + noout[i].ToString();
            }
            sw.WriteLine(dat);
            dat = it.ToString("0");
            dat = dat + "," + ttime.ToString("E");
            for(i=1;i<=nout;i++){
                dat = dat + "," + tempe[noout[i]].ToString("0.000");
            }
            sw.WriteLine(dat);
            //-----------------------------------------------------
            //時刻歴初期値入力
            //-----------------------------------------------------
            sr = new System.IO.StreamReader(ftinp, System.Text.Encoding.Default);
            dat = sr.ReadLine(); //コメント
            //*****************************************************
            //時刻歴計算（時刻歴入力ファイル：ftinpの最終行まで）
            //*****************************************************
            it = 0 ; nnn = 0;
            while(!sr.EndOfStream){
                it = it + 1;
                ttime = delta * (double)it;
                //-----------------------------------------------------
                //時刻歴入力（指定温度はその都度入力）
                //-----------------------------------------------------
                dat = sr.ReadLine() ; sbuf = dat.Split(delim); //１列目はダミー
                i = 0;
                //温度指定節点温度
                if(0<KOT){
                    for(k=1;k<=KOT;k++){
                        i = i + 1 ; Tinp[k] = double.Parse(sbuf[i]);
                    }
                }
                //熱伝達境界辺温度
                if(0<KOC){
                    for(k=1;k<=KOC;k++){
                        i = i + 1 ; Tcin2[k] = double.Parse(sbuf[i]);
                    }
                }
                //-----------------------------------------------------
                //熱流束ベクトル組立
                //-----------------------------------------------------
                //前回温度ベクトル
                for(i=1;i<=nt;i++){
                    reac[i] = 0.0;
                    for(j=1;j<=nt;j++){
                        reac[i] = reac[i] + tck1[i, j] * tempe[j];
                }}
                //温度指定境界処理
                if(0<KOT){
                    for(k=1;k<=KOT;k++){
                        for(i=1;i<=nt;i++){
                            reac[i] = reac[i] - Tinp[k] * rtck2[i, k];
                    }}
                }
                //熱伝達・発熱率ベクトル
                for(i=1;i<=nt;i++){ftvec[i] = 0.0;}
                for(ne=1;ne<=NELT;ne++){
                    //発熱率ベクトル
                    dotq = Arho[ne] * Ac[ne] * Tk[ne] * Al[ne] * Math.Exp(-Al[ne] * ttime);
                    dotq = dotq + Arho[ne] * Ac[ne] * Tk[ne] * Al[ne] * Math.Exp(-Al[ne] * (ttime - delta));
                    dotq = 0.5 * dotq;
                    FQVEC(ne, kakom, x, y, dotq, fq);
                    ASVEC(ne, 4, 1, kakom, fq, ftvec);
                    //熱伝達ベクトル
                    if(0<KOC){
                        for(k=1;k<=KOC;k++){
                            if(ne==nekc[k, 1]){
                                kchen = nekc[k, 2];
                                alpc = alphac[k];
                                tc = 0.5 * (Tcin1[k] + Tcin2[k]);
                                HTVEC(ne, kchen, kakom, x, y, alpc, tc, fc);
                                ASVEC(ne, 4, 1, kakom, fc, ftvec);
                            }
                        }
                    }
                }
                //全体ベクトル
                for(i=1;i<=nt;i++){
                    reac[i] = reac[i] + ftvec[i];
                }
                //-----------------------------------------------------
                //温度指定節点処理
                //-----------------------------------------------------
                for(i=1;i<=mm;i++){
                    ia = index[i];
                    reac[i] = reac[ia];
                }
                //-----------------------------------------------------
                //コレスキー法（前進消去・後退代入）
                //-----------------------------------------------------
                CHBAND11(mm, ib, tck2, reac);
                //-----------------------------------------------------
                //現在温度算定
                //-----------------------------------------------------
                for(i=1;i<=nt;i++){tempe[i] = 0.0 ;}
                for(i=1;i<=mm;i++){ia = index[i] ; tempe[ia] = reac[i] ;}
                for(k=1;k<=KOT;k++){tempe[nokt[k]] = Tinp[k];}
                //-----------------------------------------------------
                //記憶更新
                //-----------------------------------------------------
                if(0<KOC){
                    for(k=1;k<=KOC;k++){
                        Tcin1[k] = Tcin2[k];
                    }
                }
                //-----------------------------------------------------
                //結果出力
                //-----------------------------------------------------
                dat = it.ToString("0");
                dat = dat + "," + ttime.ToString("E");
                for(i=1;i<=nout;i++){
                    dat = dat + "," + tempe[noout[i]].ToString("0.000");
                }
                sw.WriteLine(dat);
                //-----------------------------------------------------
                //節点温度保存
                //-----------------------------------------------------
                if(0<ntout){
                    for(n=1;n<=ntout;n++){
                        if(it==nntout[n]){
                            nnn = nnn + 1;
                            for(i=1;i<=NODT;i++){
                                tempend[i, nnn] = tempe[i];
                            }
                        }
                    }
                }
            }
            sr.Close();
            //-----------------------------------------------------
            //節点温度出力
            //-----------------------------------------------------
            if(0<ntout){
                dat = "*temperature of nodes for selected times" ; sw.WriteLine(dat);
                dat = "Node,init";
                for(i=1;i<=ntout;i++){
                    dat = dat + "," + nntout[i].ToString();
                }
                sw.WriteLine(dat);
                for(i=1;i<=NODT;i++){
                    dat = i.ToString();
                    dat = dat + "," + tempe0[i].ToString("0.000");
                    for (j = 1; j <= ntout; j++)
                    {
                        dat = dat + "," + tempend[i, j].ToString("0.000");
                    }
                    sw.WriteLine(dat);
                }
            }
            //-----------------------------------------------------
            //時間計測完了
            //-----------------------------------------------------
            entime = DateTime.Now;
            tspan = entime - sttime;
            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()+" (C#)";
            sw.WriteLine(dat);
            sw.Close();
           
            label8.Text = entime.ToString();
            label10.Text = tspan.TotalSeconds.ToString("0.00sec");
            label11.Text = "nt=" + nt.ToString() + "  " + "mm=" + mm.ToString() + "  " + "ib=" + ib.ToString();
            MessageBox.Show("計算完了", "通知");
        }
        //------------------------------------------------------------------------------

        private void INPDAT(string fnameW, string strcom, int NODT, int NELT, int MATEL,
                       int KOT, int KOC, double delta, int nout,int ntout,int itmax,
                       double[] x, double[] y, double[] tempe,
                       int[] nokt, int[,] nekc, int[,] kakom,
                       double[] Ak, double[] Ac, double[] Arho, double[] Tk, double[] Al, int[] matno)
        {
            //入力確認出力
            System.IO.StreamWriter sw;
            string dat;
            int i,j,k,k1,k2,ne;

            sw = new System.IO.StreamWriter(fnameW, false, System.Text.Encoding.Default);
            dat = strcom ; sw.WriteLine(dat);
            dat = "NODT,NELT,MATEL,KOT,KOC,delta,nout,ntout,itmax" ; sw.WriteLine(dat);
            dat = NODT.ToString("0");
            dat = dat + "," + NELT.ToString("0");
            dat = dat + "," + MATEL.ToString("0");
            dat = dat + "," + KOT.ToString("0");
            dat = dat + "," + KOC.ToString("0");
            dat = dat + "," + delta.ToString("E");
            dat = dat + "," + nout.ToString("0");
            dat = dat + "," + ntout.ToString("0");
            dat = dat + "," + itmax.ToString("0");
            sw.WriteLine(dat);
            dat = "*node characteristics" ; sw.WriteLine(dat);
            dat = "node,x,y,fix-T,Tini" ; sw.WriteLine(dat);
            for(i=1;i<=NODT;i++){
                k = 0;
                if(0<KOT){
                    for(j=1;j<=KOT;j++){
                        if(i==nokt[j]){
                            k = 1 ; break;
                        }
                    }
                }
                dat = i.ToString("0");
                dat = dat + "," + x[i].ToString("E");
                dat = dat + "," + y[i].ToString("E");
                dat = dat + "," + k.ToString("0");
                dat = dat + "," + tempe[i].ToString("E");
                sw.WriteLine(dat);
            }
            dat = "*element characteristics"; sw.WriteLine(dat);
            dat = "element,node-1,node-2,node-3,node-4,k1,k2,Ak,Ac,Arho,Tk,Al,matno"; sw.WriteLine(dat);
            for(ne=1;ne<=NELT;ne++){
                k1 = 0 ; k2 = 0;
                if(0<KOC){
                    k1 = 0 ; k2 = 0;
                    for(k=1;k<=KOC;k++){
                        if(ne==nekc[k, 1]){
                            k1 = 1; //熱伝達境界を有する要素（k1=1）
                            k2 = nekc[k, 2]; //熱伝達境界の要素辺（1 or 2 or 3 or 4）
                        }
                    }   
                }
                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 + "," + kakom[ne, 4].ToString("0");
                dat = dat + "," + k1.ToString("0");
                dat = dat + "," + k2.ToString("0");
                dat = dat + "," + Ak[ne].ToString("E");
                dat = dat + "," + Ac[ne].ToString("E");
                dat = dat + "," + Arho[ne].ToString("E");
                dat = dat + "," + Tk[ne].ToString("E");
                dat = dat + "," + Al[ne].ToString("E");
                dat = dat + "," + matno[ne].ToString();
                sw.WriteLine(dat);
            }
            sw.Close();
        }
        //------------------------------------------------------------------------------
        private void NMAT(int ne, int[,] kakom, double[] x, double[] y,
                     double[,] dnk, double[,] dnc, out double detJ, double a, double b)
        {
            int i,j,k,l;
            double J11,J12,J21,J22;
            double dn1a,dn2a,dn3a,dn4a,dn1b,dn2b,dn3b,dn4b;
            double dn1x,dn2x,dn3x,dn4x,dn1y,dn2y,dn3y,dn4y;
            double nm1,nm2,nm3,nm4;

            i = kakom[ne, 1];
            j = kakom[ne, 2];
            k = kakom[ne, 3];
            l = kakom[ne, 4];

            //[dN/da dN/db]
            dn1a = -0.25 * (1.0 - b);
            dn2a =  0.25 * (1.0 - b);
            dn3a =  0.25 * (1.0 + b);
            dn4a = -0.25 * (1.0 + b);
            dn1b = -0.25 * (1.0 - a);
            dn2b = -0.25 * (1.0 + a);
            dn3b =  0.25 * (1.0 + a);
            dn4b =  0.25 * (1.0 - a);
            //Jacobiマトリックスとdet[J]
            J11 = dn1a * x[i] + dn2a * x[j] + dn3a * x[k] + dn4a * x[l];
            J12 = dn1a * y[i] + dn2a * y[j] + dn3a * y[k] + dn4a * y[l];
            J21 = dn1b * x[i] + dn2b * x[j] + dn3b * x[k] + dn4b * x[l];
            J22 = dn1b * y[i] + dn2b * y[j] + dn3b * y[k] + dn4b * y[l];
            detJ = J11 * J22 - J12 * J21;
            //[dN/dx dN/dy]
            dn1x = ( J22 * dn1a - J12 * dn1b) / detJ;
            dn2x = ( J22 * dn2a - J12 * dn2b) / detJ;
            dn3x = ( J22 * dn3a - J12 * dn3b) / detJ;
            dn4x = ( J22 * dn4a - J12 * dn4b) / detJ;
            dn1y = (-J21 * dn1a + J11 * dn1b) / detJ;
            dn2y = (-J21 * dn2a + J11 * dn2b) / detJ;
            dn3y = (-J21 * dn3a + J11 * dn3b) / detJ;
            dn4y = (-J21 * dn4a + J11 * dn4b) / detJ;
            //[N]行列要素
            nm1 = 0.25 * (1.0 - a) * (1.0 - b);
            nm2 = 0.25 * (1.0 + a) * (1.0 - b);
            nm3 = 0.25 * (1.0 + a) * (1.0 + b);
            nm4 = 0.25 * (1.0 - a) * (1.0 + b);

            dnk[1, 1] = dn1x * dn1x + dn1y * dn1y;
            dnk[1, 2] = dn1x * dn2x + dn1y * dn2y;
            dnk[1, 3] = dn1x * dn3x + dn1y * dn3y;
            dnk[1, 4] = dn1x * dn4x + dn1y * dn4y;
            dnk[2, 1] = dn2x * dn1x + dn2y * dn1y;
            dnk[2, 2] = dn2x * dn2x + dn2y * dn2y;
            dnk[2, 3] = dn2x * dn3x + dn2y * dn3y;
            dnk[2, 4] = dn2x * dn4x + dn2y * dn4y;
            dnk[3, 1] = dn3x * dn1x + dn3y * dn1y;
            dnk[3, 2] = dn3x * dn2x + dn3y * dn2y;
            dnk[3, 3] = dn3x * dn3x + dn3y * dn3y;
            dnk[3, 4] = dn3x * dn4x + dn3y * dn4y;
            dnk[4, 1] = dn4x * dn1x + dn4y * dn1y;
            dnk[4, 2] = dn4x * dn2x + dn4y * dn2y;
            dnk[4, 3] = dn4x * dn3x + dn4y * dn3y;
            dnk[4, 4] = dn4x * dn4x + dn4y * dn4y;

            dnc[1, 1] = nm1 * nm1;
            dnc[1, 2] = nm1 * nm2;
            dnc[1, 3] = nm1 * nm3;
            dnc[1, 4] = nm1 * nm4;
            dnc[2, 1] = nm2 * nm1;
            dnc[2, 2] = nm2 * nm2;
            dnc[2, 3] = nm2 * nm3;
            dnc[2, 4] = nm2 * nm4;
            dnc[3, 1] = nm3 * nm1;
            dnc[3, 2] = nm3 * nm2;
            dnc[3, 3] = nm3 * nm3;
            dnc[3, 4] = nm3 * nm4;
            dnc[4, 1] = nm4 * nm1;
            dnc[4, 2] = nm4 * nm2;
            dnc[4, 3] = nm4 * nm3;
            dnc[4, 4] = nm4 * nm4;
        }
        //------------------------------------------------------------------------------
        private void IPOINT(int kk, out double a, out double b)
        {
           //面積積分におけるGauss積分点
            a=0.0;b=0.0;
            switch(kk){
                case 1 : a = -0.5773502692 ; b = -0.5773502692;break;
                case 2 : a =  0.5773502692 ; b = -0.5773502692;break;
                case 3 : a =  0.5773502692 ; b =  0.5773502692;break;
                case 4 : a = -0.5773502692 ; b =  0.5773502692;break;
            }
        }
        //------------------------------------------------------------------------------
        private void STIFF(int ne, int[,] kakom, double[] x, double[] y,
                      double[] Ak, double[] Ac, double[] Arho, double[,] ck1, double[,] ck2, double delta)
        {
            //要素剛性行列
            int kk,i,j;
            double[,] dnk=new double[5, 5];
            double[,] dnc=new double[5, 5];
            double[,] ww1=new double[5, 5];
            double[,] ww2=new double[5, 5];
            double detJ,a,b;

            for(i=1;i<=4;i++){
                for(j=1;j<=4;j++){
                    ww1[i, j] = 0.0;
                    ww2[i, j] = 0.0;
                    ck1[i, j] = 0.0;
                    ck2[i, j] = 0.0;
            }}
            for(kk=1;kk<=4;kk++){
                IPOINT(kk, out a, out b);
                NMAT(ne, kakom, x, y, dnk, dnc, out detJ, a, b);
                for(i=1;i<=4;i++){
                    for(j=1;j<=4;j++){
                        ww1[i, j] = ww1[i, j] + Ak[ne] * dnk[i, j] * detJ;
                        ww2[i, j] = ww2[i, j] + Arho[ne] * Ac[ne] * dnc[i, j] * detJ;
                }}
            }
            for(i=1;i<=4;i++){
                for(j=1;j<=4;j++){
                    ck1[i, j] = -0.5 * ww1[i, j] + 1.0 / delta * ww2[i, j];
                    ck2[i, j] = +0.5 * ww1[i, j] + 1.0 / delta * ww2[i, j];
            }}
        }
        //------------------------------------------------------------------------------
        private void LPOINT(int kk, int kchen, out double a, out double b)
        {
            //要素辺上積分におけるGauss積分点
            a=0.0;b=0.0;
            switch(kchen){
                case 1:
                    b = -1.0;
                    if( kk == 1 ) a = -0.5773502692;
                    if( kk == 2 ) a = 0.5773502692;
                    break;
                case 2:
                    a = 1.0;
                    if( kk == 1 ) b = -0.5773502692;
                    if( kk == 2 ) b = 0.5773502692;
                    break;
                case 3:
                    b = 1.0;
                    if( kk == 1 ) a = -0.5773502692;
                    if( kk == 2 ) a = 0.5773502692;
                    break;
                case 4:
                    a = -1.0;
                    if( kk == 1 ) b = -0.5773502692;
                    if( kk == 2 ) b = 0.5773502692;
                    break;
            }
        }
        //------------------------------------------------------------------------------
        private void HTMAT(int ne, int kchen, int[,] kakom, double[] x, double[] y,
                          double alpc, double[,] ht)
        {
            //要素熱伝達マトリックス
            int i,j,kk;
            double a,b,s;
            double[] vn=new double[5];

            i = 0; j = 0;
            switch(kchen){
                case 1 : i = kakom[ne, 1] ; j = kakom[ne, 2];break;
                case 2 : i = kakom[ne, 2] ; j = kakom[ne, 3];break;
                case 3 : i = kakom[ne, 3] ; j = kakom[ne, 4];break;
                case 4 : i = kakom[ne, 4] ; j = kakom[ne, 1];break;
            }
            s = Math.Sqrt((x[j]-x[i])*(x[j]-x[i]) + (y[j]-y[i])*(y[j]-y[i]));

            for(i=1;i<=4;i++){
                for(j=1;j<=4;j++){
                    ht[i, j] = 0.0;
            }}
            for(kk=1;kk<=2;kk++){
                LPOINT(kk, kchen, out a, out b);
                //[N]行列要素
                vn[1] = 0.25 * (1.0 - a) * (1.0 - b);
                vn[2] = 0.25 * (1.0 + a) * (1.0 - b);
                vn[3] = 0.25 * (1.0 + a) * (1.0 + b);
                vn[4] = 0.25 * (1.0 - a) * (1.0 + b);
                for(i=1;i<=4;i++){
                    for(j=1;j<=4;j++){
                        ht[i, j] = ht[i, j] + 0.5 * s * alpc * vn[i] * vn[j];
                }}
            }
        }
        //------------------------------------------------------------------------------
        private void FQVEC(int ne, int[,] kakom, double[] x, double[] y, double dotq, double[] fq)
        {
            //発熱率ベクトル
            int i,j,k,l,kk;
            double a,b,detJ;
            double[] vn=new double[5];
            double J11,J12,J21,J22;
            double dn1a,dn2a,dn3a,dn4a,dn1b,dn2b,dn3b,dn4b;

            for(i=1;i<=4;i++){fq[i] = 0.0;}
            for(kk=1;kk<=4;kk++){
                IPOINT(kk, out a, out b);
                i = kakom[ne, 1];
                j = kakom[ne, 2];
                k = kakom[ne, 3];
                l = kakom[ne, 4];
                //[dN/da dN/db]
                dn1a = -0.25 * (1.0 - b);
                dn2a =  0.25 * (1.0 - b);
                dn3a =  0.25 * (1.0 + b);
                dn4a = -0.25 * (1.0 + b);
                dn1b = -0.25 * (1.0 - a);
                dn2b = -0.25 * (1.0 + a);
                dn3b =  0.25 * (1.0 + a);
                dn4b =  0.25 * (1.0 - a);
                //Jacobiマトリックスとdet[J]
                J11 = dn1a * x[i] + dn2a * x[j] + dn3a * x[k] + dn4a * x[l];
                J12 = dn1a * y[i] + dn2a * y[j] + dn3a * y[k] + dn4a * y[l];
                J21 = dn1b * x[i] + dn2b * x[j] + dn3b * x[k] + dn4b * x[l];
                J22 = dn1b * y[i] + dn2b * y[j] + dn3b * y[k] + dn4b * y[l];
                detJ = J11 * J22 - J12 * J21;
                //[N]行列要素
                vn[1] = 0.25 * (1.0 - a) * (1.0 - b);
                vn[2] = 0.25 * (1.0 + a) * (1.0 - b);
                vn[3] = 0.25 * (1.0 + a) * (1.0 + b);
                vn[4] = 0.25 * (1.0 - a) * (1.0 + b);
                for(i=1;i<=4;i++){
                    fq[i] = fq[i] + dotq * vn[i] * detJ;
                }
            }
        }
        //------------------------------------------------------------------------------
        private void HTVEC(int ne, int kchen, int[,] kakom, double[] x, double[] y,
                          double alpc, double tc, double[] fc)
        {
            //熱伝達ベクトル
            int i,j,kk;
            double[] vn=new double[5];
            double a,b,s;

            i = 0; j = 0;
            switch(kchen){
                case 1 : i = kakom[ne, 1] ; j = kakom[ne, 2];break;
                case 2 : i = kakom[ne, 2] ; j = kakom[ne, 3];break;
                case 3 : i = kakom[ne, 3] ; j = kakom[ne, 4];break;
                case 4 : i = kakom[ne, 4] ; j = kakom[ne, 1];break;
            }
            s = Math.Sqrt((x[j] - x[i])*(x[j] - x[i]) + (y[j] - y[i])*(y[j] - y[i]));

            for(i=1;i<=4;i++){fc[i] = 0.0;}
            for(kk=1;kk<=2;kk++){
                LPOINT(kk, kchen, out a, out b);
                //[N]行列要素
                vn[1] = 0.25 * (1.0 - a) * (1.0 - b);
                vn[2] = 0.25 * (1.0 + a) * (1.0 - b);
                vn[3] = 0.25 * (1.0 + a) * (1.0 + b);
                vn[4] = 0.25 * (1.0 - a) * (1.0 + b);
                for(i=1;i<=4;i++){
                    fc[i] = fc[i] + 0.5 * s * alpc * tc * vn[i];
                }
            }
        }
        //------------------------------------------------------------------------------
        private double CALtempeNEL(int ne, int[,] kakom, double[] tempe)
        {
            //要素平均温度算定
            int i,j,k,l,kk;
            double nm1,nm2,nm3,nm4,a,b,tem;

            i = kakom[ne, 1];
            j = kakom[ne, 2];
            k = kakom[ne, 3];
            l = kakom[ne, 4];
            tem = 0.0;
            for(kk=1;kk<=4;kk++){
                IPOINT(kk, out a, out b);
                //[N]行列要素
                nm1 = 0.25 * (1.0 - a) * (1.0 - b);
                nm2 = 0.25 * (1.0 + a) * (1.0 - b);
                nm3 = 0.25 * (1.0 + a) * (1.0 + b);
                nm4 = 0.25 * (1.0 - a) * (1.0 + b);
                tem = tem + nm1 * tempe[i] + nm2 * tempe[j] + nm3 * tempe[k] + nm4 * tempe[l];
            }
            return 0.25 * tem;
        }
        //------------------------------------------------------------------------------
        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 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();
        }
        //------------------------------------------------------------------------------

    }
}