Option Explicit On
Option Strict On
Public Class Form1

    Private Sub Form1_Load(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles MyBase.Load
        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"
    End Sub


    Private Sub ToolStripButton1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles ToolStripButton1.Click
        Call MAIN()
    End Sub

    Private Sub MAIN()
        '平面骨組構造解析
        Dim i As Integer : Dim j As Integer
        Dim ne As Integer : Dim n As Integer
        Dim ia As Integer : Dim ja As Integer
        Dim sr As System.IO.StreamReader
        Dim sw As System.IO.StreamWriter
        Dim dat As String
        Dim sbuf() As String
        Dim delim() As Char = {","c}
        Dim fnameR As String = ""
        Dim fnameW As String = ""
        Dim strcom As String

        Dim NODT As Integer       '節点総数
        Dim NELT As Integer       '要素総数
        Dim MATEL As Integer      '材料種別数
        Dim KOX As Integer        'x方向変位既知節点数
        Dim KOY As Integer        'y方向変位既知節点数
        Dim KOZ As Integer        '回転変位既知節点数
        Dim NF As Integer         '載荷節点数
        Dim nt As Integer         '全自由度
        Dim mm As Integer         '変位既知節点処理語の自由度
        Dim kakom(,) As Integer   '要素構成接点番号
        Dim Em() As Double        '要素弾性係数
        Dim AA() As Double        '要素断面積
        Dim AI() As Double        '要素断面二次モーメント
        Dim gamma() As Double     '要素密度
        Dim gkh() As Double       '要素水平加速度
        Dim gkv() As Double       '要素鉛直加速度
        Dim alpha() As Double     '要素熱膨張係数
        Dim matno() As Integer    '材料種別No
        Dim x() As Double         '節点x座標
        Dim y() As Double         '節点y座標
        Dim deltaT() As Double    '節点温度変化量
        Dim nokx() As Integer     'x方向変位既知節点番号
        Dim noky() As Integer     'y方向変位既知節点番号
        Dim nokz() As Integer     '回転変位既知節点番号
        Dim f() As Double         '全体節点外力ベクトル（入力値）
        Dim fmass() As Double     '全体慣性力ベクトル（入力値）
        Dim ftemp() As Double    '全体温度荷重相当力ベクトル（入力値）
        Dim ftvec() As Double     '全体節点力ベクトル
        Dim fdis() As Double      '既知節点変位による荷重成分格納
        Dim rdis() As Double      '既知節点変位量
        Dim index() As Integer    '変位既知節点処理のためのインデックス
        Dim sm(6, 6) As Double    '要素剛性マトリックス
        Dim tsm(,) As Double      '全体剛性マトリックス
        Dim dis() As Double       '全体節点変位
        Dim reac() As Double      '要素節点反力
        Dim hen(6, 6) As Double   '座標変換行列
        Dim fevec(6) As Double    '要素節点力ベクトル
        Dim wm(,) As Double       '材料種別作業用
        Dim sttime As Date        '計算開始時刻
        Dim entime As Date        '計算終了時刻
        Dim tspan As TimeSpan     '計算時間
        Dim ib As Integer         'バンド幅
        Dim fact0 As Double = 1.0E-30 '0判定
        Dim ggg As Double = 9.8   '重力加速度
        '重力加速度は9.8に固定．慣性力ベクトルを求める中でgで除してgを乗じるため値は任意でかまわない

        '-----------------------------------------------------
        '入出力データファイル指定
        '-----------------------------------------------------
        OpenFileDialog1.InitialDirectory = My.Computer.FileSystem.CurrentDirectory
        If OpenFileDialog1.ShowDialog() = Windows.Forms.DialogResult.OK Then fnameR = OpenFileDialog1.FileName
        If SaveFileDialog1.ShowDialog() = Windows.Forms.DialogResult.OK Then fnameW = SaveFileDialog1.FileName
        Label2.Text = fnameR
        Label4.Text = fnameW
        '-----------------------------------------------------
        '時間計測開始
        '-----------------------------------------------------
        sttime = DateTime.Now
        Label6.Text = sttime.ToString
        '*****************************************************
        'データ入力
        '*****************************************************
        sr = New System.IO.StreamReader(fnameR, System.Text.Encoding.Default)
        '--------------------------------------------------------------------------------------
        '書出用コメント入力
        '--------------------------------------------------------------------------------------
        dat = sr.ReadLine : sbuf = dat.Split(delim) : strcom = sbuf(0)
        '--------------------------------------------------------------------------------------
        '節点数，要素数，材料種別数,X方向変位既知節点数，Y方向変位既知節点数，回転変位既知節点数，載荷点数
        '--------------------------------------------------------------------------------------
        dat = sr.ReadLine : sbuf = dat.Split(delim)
        NODT = CInt(sbuf(0))
        NELT = CInt(sbuf(1))
        MATEL = CInt(sbuf(2))
        KOX = CInt(sbuf(3))
        KOY = CInt(sbuf(4))
        KOZ = CInt(sbuf(5))
        NF = CInt(sbuf(6))
        nt = NODT * 3
        '--------------------------------------------------------------------------------------
        '配列寸法宣言
        '--------------------------------------------------------------------------------------
        ReDim kakom(NELT, 2)
        ReDim Em(NELT)
        ReDim AA(NELT)
        ReDim AI(NELT)
        ReDim gamma(NELT)
        ReDim gkh(NELT)
        ReDim gkv(NELT)
        ReDim alpha(NELT)
        ReDim matno(NELT)
        ReDim wm(MATEL, 7)
        ReDim x(NODT)
        ReDim y(NODT)
        ReDim deltaT(NODT)
        ReDim nokx(KOX)
        ReDim noky(KOY)
        ReDim nokz(KOZ)
        ReDim index(nt)
        ReDim tsm(nt, nt)
        ReDim f(nt)
        ReDim fmass(nt)
        ReDim ftemp(nt)
        ReDim ftvec(nt)
        ReDim fdis(nt)
        ReDim rdis(nt)
        ReDim dis(nt)
        ReDim reac(nt)
        '------------------------------------------------------------------------------------------
        '材料特性（Em,AA,AI）
        '------------------------------------------------------------------------------------------
        For i = 1 To MATEL
            dat = sr.ReadLine : sbuf = dat.Split(delim)
            For j = 1 To 7
                wm(i, j) = CDbl(sbuf(j - 1))
            Next j
        Next i
        '------------------------------------------------------------------------------------------
        '要素構成節点と材料種別（node-1,node-2,matno）
        '------------------------------------------------------------------------------------------
        For ne = 1 To NELT
            dat = sr.ReadLine : sbuf = dat.Split(delim)
            kakom(ne, 1) = CInt(sbuf(0))
            kakom(ne, 2) = CInt(sbuf(1))
            matno(ne) = CInt(sbuf(2))
            For i = 1 To MATEL
                If i = matno(ne) Then
                    Em(ne) = wm(i, 1)
                    AA(ne) = wm(i, 2)
                    AI(ne) = wm(i, 3)
                    gamma(ne) = wm(i, 4)
                    gkh(ne) = wm(i, 5)
                    gkv(ne) = wm(i, 6)
                    alpha(ne) = wm(i, 7)
                End If
            Next i
        Next ne
        '--------------------------------------------------------------------------------------
        '節点座標（x,y）
        '--------------------------------------------------------------------------------------
        For i = 1 To NODT
            dat = sr.ReadLine : sbuf = dat.Split(delim)
            x(i) = CDbl(sbuf(0))
            y(i) = CDbl(sbuf(1))
            deltaT(i) = CDbl(sbuf(2))
        Next i
        '--------------------------------------------------------------------------------------
        '変位既知節点番号，既知節点変位量
        '--------------------------------------------------------------------------------------
        For i = 1 To nt
            rdis(i) = 0.0
        Next i
        If 0 < KOX Then
            For i = 1 To KOX
                dat = sr.ReadLine : sbuf = dat.Split(delim)
                nokx(i) = CInt(sbuf(0))
                rdis(3 * nokx(i) - 2) = CDbl(sbuf(1))
            Next i
        End If
        If 0 < KOY Then
            For i = 1 To KOY
                dat = sr.ReadLine : sbuf = dat.Split(delim)
                noky(i) = CInt(sbuf(0))
                rdis(3 * noky(i) - 1) = CDbl(sbuf(1))
            Next i
        End If
        If 0 < KOZ Then
            For i = 1 To KOZ
                dat = sr.ReadLine : sbuf = dat.Split(delim)
                nokz(i) = CInt(sbuf(0))
                rdis(3 * nokz(i)) = CDbl(sbuf(1))
            Next i
        End If
        '--------------------------------------------------------------------------------------
        '節点番号,x方向荷重,y方向荷重,モーメント荷重
        '--------------------------------------------------------------------------------------
        For i = 1 To nt : f(i) = 0.0 : Next i
        If 0 < NF Then
            For i = 1 To NF
                dat = sr.ReadLine : sbuf = dat.Split(delim)
                n = CInt(sbuf(0))
                f(3 * n - 2) = CDbl(sbuf(1))
                f(3 * n - 1) = CDbl(sbuf(2))
                f(3 * n) = CDbl(sbuf(3))
            Next i
        End If
        sr.Close()
        '*****************************************************
        '入力確認出力
        '*****************************************************
        Call INPDAT(fnameW, strcom, NODT, NELT, MATEL, KOX, KOY, KOZ, NF, x, y, deltaT, f, _
                    rdis, nokx, noky, nokz, kakom, Em, AA, AI, gamma, gkh, gkv, alpha, matno)
        '*****************************************************
        'メイン処理計算
        '*****************************************************
        '-----------------------------------------------------
        '全体剛性行列・全体荷重ベクトルクリア
        '-----------------------------------------------------
        For i = 1 To nt
            ftvec(i) = 0.0
            fmass(i) = 0.0
            ftemp(i) = 0.0
            dis(i) = 0.0
            reac(i) = 0.0
            For j = 1 To nt
                tsm(i, j) = 0.0
        Next j, i
        '-----------------------------------------------------
        '要素剛性行列組立
        '-----------------------------------------------------

        For ne = 1 To NELT
            '剛性行列組立
            Call SM_FRM(ne, kakom, x, y, Em, AA, AI, sm)     '要素剛性行列
            Call GLTMX(ne, kakom, x, y, sm, hen)            '座標変換行列と[ak]の座標変換
            Call ASMAT(ne, 2, 3, kakom, sm, tsm)           '全体行列組立：[ak]→[TAK]

            '慣性力ベクトル組立
            Call MM_FRM(ne, kakom, x, y, ggg, AA, gamma, sm) '要素質量行列
            Call GLTMX(ne, kakom, x, y, sm, hen)            '座標変換行列と[am]の座標変換
            Call MV_FRM(ne, sm, ggg, gkh, gkv, fevec)
            Call ASVEC(ne, 2, 3, kakom, fevec, fmass)

            '温度荷重相当ベクトル組立
            Call TM_FRM(ne, kakom, deltaT, Em, AA, alpha, fevec, hen)
            Call ASVEC(ne, 2, 3, kakom, fevec, ftemp)
        Next ne
        '-----------------------------------------------------
        '既知節点変位による荷重成分ベクトル作成
        '-----------------------------------------------------
        For i = 1 To nt
            fdis(i) = 0.0
            For j = 1 To nt
                fdis(i) = fdis(i) + rdis(j) * tsm(i, j)
        Next j, i
        '-----------------------------------------------------
        '荷重ベクトルセット
        '-----------------------------------------------------
        For i = 1 To nt
            ftvec(i) = f(i) + fmass(i) + ftemp(i) - fdis(i)
        Next i
        '-----------------------------------------------------
        '変位既知節点自由度記憶
        '-----------------------------------------------------
        For i = 1 To nt
            index(i) = i
            reac(i) = ftvec(i)
        Next i
        If 0 < KOX Then
            For i = 1 To KOX : n = 3 * nokx(i) - 2 : index(n) = 0 : Next i
        End If
        If 0 < KOY Then
            For i = 1 To KOY : n = 3 * noky(i) - 1 : index(n) = 0 : Next i
        End If
        If 0 < KOZ Then
            For i = 1 To KOZ : n = 3 * nokz(i) : index(n) = 0 : Next i
        End If
        mm = 0
        For i = 1 To nt
            If index(i) > 0 Then
                mm = mm + 1 : index(mm) = i
            End If
        Next i
        '-----------------------------------------------------
        '変位既知節点処理（連立一次方程式縮小）
        '-----------------------------------------------------
        For i = 1 To mm
            ia = index(i)
            reac(i) = reac(ia)
            For j = 1 To mm
                ja = index(j)
                tsm(i, j) = tsm(ia, ja)
        Next j, i
        '-----------------------------------------------------
        'バンド幅計算とフルマトリックスのバンド化記憶
        '-----------------------------------------------------
        ib = IBAND(mm, tsm, fact0)
        Call BAND(mm, ib, tsm)
        '-----------------------------------------------------
        '対角項確認と異常終了通知
        '-----------------------------------------------------
        For i = 1 To mm
            If Math.Abs(tsm(i, 1)) < fact0 Then Call IERROR2(fnameW)
        Next i
        '-----------------------------------------------------
        '連立一次方程式（コレスキー法）
        '-----------------------------------------------------
        Call CHBAND10(mm, ib, tsm)       '三角行列
        Call CHBAND11(mm, ib, tsm, reac) '前進消去・後退代入
        '-----------------------------------------------------
        '変位計算
        '-----------------------------------------------------
        For i = 1 To nt : dis(i) = 0.0 : Next i
        For i = 1 To mm
            ia = index(i) : dis(ia) = reac(i) '連立方程式の解としての変位セット
        Next i
        For i = 1 To nt
            reac(i) = 0.0
            dis(i) = dis(i) + rdis(i) '既知節点変位の足し込み
        Next i
        '-----------------------------------------------------
        '断面力計算
        '-----------------------------------------------------
        For ne = 1 To NELT
            Call SM_FRM(ne, kakom, x, y, Em, AA, AI, sm)     '要素剛性行列
            Call GLTMX(ne, kakom, x, y, sm, hen)            '座標変換行列と[ak]の座標変換
            Call CALST(ne, kakom, dis, sm, hen, reac, tsm, Em, AA, alpha, deltaT)
        Next ne
        '-----------------------------------------------------
        '時間計測完了
        '-----------------------------------------------------
        entime = DateTime.Now
        tspan = entime - sttime
        Label8.Text = entime.ToString
        Label10.Text = tspan.TotalSeconds.ToString("0.00sec")
        Label11.Text = "NODT=" & NODT.ToString() & "  nt=" & nt.ToString() & "  mm=" & mm.ToString() & "  ib=" & ib.ToString()
        '*****************************************************
        '結果出力
        '*****************************************************
        Call OUTRESUL(fnameW, NODT, NELT, reac, ftvec, x, y, dis, tsm, matno)

        sw = New System.IO.StreamWriter(fnameW, True, System.Text.Encoding.Default)
        dat = "NODT=" & NODT.ToString() & "  nt=" & nt.ToString() & "  mm=" & mm.ToString() & "  ib=" & ib.ToString()
        sw.WriteLine(dat)
        dat = "Calculation time=" & tspan.TotalSeconds.ToString("0.00sec")
        sw.WriteLine(dat)
        dat = entime.ToString
        sw.WriteLine(dat)
        sw.Close()

        MessageBox.Show("計算完了", "通知")

    End Sub


    Private Sub INPDAT(ByRef fnameW As String, ByVal strcom As String, _
                   ByVal NODT As Integer, ByVal NELT As Integer, ByVal MATEL As Integer, _
                   ByVal KOX As Integer, ByVal KOY As Integer, ByVal KOZ As Integer, _
                   ByVal NF As Integer, ByVal x() As Double, ByRef y() As Double, ByRef deltaT() As Double, _
                   ByRef f() As Double, ByRef rdis() As Double, _
                   ByRef nokx() As Integer, ByRef noky() As Integer, ByRef nokz() As Integer, _
                   ByRef kakom(,) As Integer, ByRef Em() As Double, ByRef AA() As Double, ByRef AI() As Double, _
                   ByRef gamma() As Double, ByRef gkh() As Double, ByRef gkv() As Double, ByRef alpha() As Double, _
                   ByRef matno() As Integer)

        '入力確認出力
        Dim sw As System.IO.StreamWriter
        Dim dat As String
        Dim i As Integer
        Dim k As Integer
        Dim l As Integer
        Dim m As Integer
        Dim ne As Integer
        sw = New System.IO.StreamWriter(fnameW, False, System.Text.Encoding.Default)
        dat = strcom : sw.WriteLine(dat)
        dat = "NODT,NELT,MATEL,KOX,KOY,KOZ,NF" : sw.WriteLine(dat)
        dat = NODT.ToString("0")
        dat = dat & "," & NELT.ToString("0")
        dat = dat & "," & MATEL.ToString("0")
        dat = dat & "," & KOX.ToString("0")
        dat = dat & "," & KOY.ToString("0")
        dat = dat & "," & KOZ.ToString("0")
        dat = dat & "," & NF.ToString("0")
        sw.WriteLine(dat)
        dat = "*node characteristics" : sw.WriteLine(dat)
        dat = "node,x,y,fx,fy,fz,fix-x,fix-y,fix-z,rdis-x,rdis-y,rdis-z,deltaT" : sw.WriteLine(dat)
        For i = 1 To NODT
            k = 0 : l = 0 : m = 0
            If 0 < KOX Then
                For j = 1 To KOX
                    If i = nokx(j) Then k = 1
                Next j
            End If
            If 0 < KOY Then
                For j = 1 To KOY
                    If i = noky(j) Then l = 1
                Next j
            End If
            If 0 < KOZ Then
                For j = 1 To KOZ
                    If i = nokz(j) Then m = 1
                Next j
            End If
            dat = i.ToString("0")
            dat = dat & "," & x(i).ToString("E")
            dat = dat & "," & y(i).ToString("E")
            dat = dat & "," & f(3 * i - 2).ToString("E")
            dat = dat & "," & f(3 * i - 1).ToString("E")
            dat = dat & "," & f(3 * i).ToString("E")
            dat = dat & "," & k.ToString("0")
            dat = dat & "," & l.ToString("0")
            dat = dat & "," & m.ToString("0")
            dat = dat & "," & rdis(3 * i - 2).ToString("E")
            dat = dat & "," & rdis(3 * i - 1).ToString("E")
            dat = dat & "," & rdis(3 * i).ToString("E")
            dat = dat & "," & deltaT(i).ToString("E")
            sw.WriteLine(dat)
        Next i
        dat = "*element characteristics" : sw.WriteLine(dat)
        dat = "element,node-1,node-2,E,A,I,gamma,gkh,gkv,alpha,matno" : sw.WriteLine(dat)
        For ne = 1 To NELT
            dat = ne.ToString("0")
            dat = dat & "," & kakom(ne, 1).ToString("0")
            dat = dat & "," & kakom(ne, 2).ToString("0")
            dat = dat & "," & Em(ne).ToString("E")
            dat = dat & "," & AA(ne).ToString("E")
            dat = dat & "," & AI(ne).ToString("E")
            dat = dat & "," & gamma(ne).ToString("E")
            dat = dat & "," & gkh(ne).ToString("E")
            dat = dat & "," & gkv(ne).ToString("E")
            dat = dat & "," & alpha(ne).ToString("E")
            dat = dat & "," & matno(ne).ToString("0")
            sw.WriteLine(dat)
        Next ne
        sw.Close()

    End Sub


    Private Sub OUTRESUL(ByVal fnameW As String, ByVal NODT As Integer, ByVal NELT As Integer, _
                     ByRef reac() As Double, ByRef ftvec() As Double, _
                     ByRef x() As Double, ByRef y() As Double, ByRef dis() As Double, _
                     ByRef tsm(,) As Double, ByRef matno() As Integer)
        '結果出力（Appendモード）
        Dim i As Integer
        Dim ne As Integer
        Dim sw As System.IO.StreamWriter
        Dim dat As String

        sw = New System.IO.StreamWriter(fnameW, True, System.Text.Encoding.Default)
        dat = "*displacements and forces" : sw.WriteLine(dat)
        dat = "node,coord-x,coord-y,dir-x,dir-y,dir-z,reac-x,reac-y,reac-z,ftvec-x,ftvec-y,ftvec-z" : sw.WriteLine(dat)
        For i = 1 To NODT
            dat = i.ToString("0")
            dat = dat & "," & x(i).ToString("E")
            dat = dat & "," & y(i).ToString("E")
            dat = dat & "," & dis(3 * i - 2).ToString("E")
            dat = dat & "," & dis(3 * i - 1).ToString("E")
            dat = dat & "," & dis(3 * i).ToString("E")
            dat = dat & "," & reac(3 * i - 2).ToString("E")
            dat = dat & "," & reac(3 * i - 1).ToString("E")
            dat = dat & "," & reac(3 * i).ToString("E")
            dat = dat & "," & ftvec(3 * i - 2).ToString("E")
            dat = dat & "," & ftvec(3 * i - 1).ToString("E")
            dat = dat & "," & ftvec(3 * i).ToString("E")
            sw.WriteLine(dat)
        Next i
        dat = "*stress resultants" : sw.WriteLine(dat)
        dat = "element,Ni,Si,Mi,Nj,Sj,Mj" : sw.WriteLine(dat)
        For ne = 1 To NELT
            dat = ne.ToString("0")
            For i = 1 To 6
                dat = dat & "," & tsm(ne, i).ToString("E")
            Next i
            sw.WriteLine(dat)
        Next ne
        sw.Close()

    End Sub


    Private Sub SM_FRM(ByVal ne As Integer, ByRef kakom(,) As Integer, ByRef x() As Double, ByRef y() As Double, _
                    ByRef Em() As Double, ByRef AA() As Double, ByRef AI() As Double, ByRef ak(,) As Double)
        '要素剛性行列
        Dim i As Integer
        Dim j As Integer
        Dim AL As Double
        Dim EA As Double
        Dim EI As Double

        '行列初期化
        For i = 1 To 6
            For j = 1 To 6
                ak(i, j) = 0.0
        Next j, i

        '要素長さ
        i = kakom(ne, 1)
        j = kakom(ne, 2)
        AL = Math.Sqrt((x(j) - x(i)) ^ 2 + (y(j) - y(i)) ^ 2)

        '要素剛性行列
        EA = Em(ne) * AA(ne)
        EI = Em(ne) * AI(ne)
        ak(1, 1) = EA / AL
        ak(2, 1) = 0.0
        ak(3, 1) = 0.0
        ak(4, 1) = -EA / AL
        ak(5, 1) = 0.0
        ak(6, 1) = 0.0
        ak(2, 2) = 12.0 * EI / AL / AL / AL
        ak(3, 2) = 6.0 * EI / AL / AL
        ak(4, 2) = 0.0
        ak(5, 2) = -12.0 * EI / AL / AL / AL
        ak(6, 2) = 6.0 * EI / AL / AL
        ak(3, 3) = 4.0 * EI / AL
        ak(4, 3) = 0.0
        ak(5, 3) = -6.0 * EI / AL / AL
        ak(6, 3) = 2.0 * EI / AL
        ak(4, 4) = EA / AL
        ak(5, 4) = 0.0
        ak(6, 4) = 0.0
        ak(5, 5) = 12.0 * EI / AL / AL / AL
        ak(6, 5) = -6.0 * EI / AL / AL
        ak(6, 6) = 4.0 * EI / AL
        '対称項セット
        For i = 1 To 5
            For j = i + 1 To 6
                ak(i, j) = ak(j, i)
            Next j
        Next i

    End Sub


    Private Sub MM_FRM(ByVal ne As Integer, ByRef kakom(,) As Integer, ByRef x() As Double, ByRef y() As Double, _
                        ByVal ggg As Double, ByRef AA() As Double, ByRef gamma() As Double, ByRef am(,) As Double)
        '要素質量行列
        Dim i As Integer
        Dim j As Integer
        Dim AL As Double

        '行列初期化
        For i = 1 To 6
            For j = 1 To 6
                am(i, j) = 0.0
        Next j, i

        '要素長さ
        i = kakom(ne, 1)
        j = kakom(ne, 2)
        AL = Math.Sqrt((x(j) - x(i)) ^ 2 + (y(j) - y(i)) ^ 2)

        am(1, 1) = 1.0 / 3.0
        am(2, 1) = 0.0
        am(3, 1) = 0.0
        am(4, 1) = 1.0 / 6.0
        am(5, 1) = 0.0
        am(6, 1) = 0.0
        am(2, 2) = 13.0 / 35.0
        am(3, 2) = 11.0 * AL / 210.0
        am(4, 2) = 0.0
        am(5, 2) = 9.0 / 70.0
        am(6, 2) = -13.0 * AL / 420.0
        am(3, 3) = AL * AL / 105.0
        am(4, 3) = 0.0
        am(5, 3) = 13.0 * AL / 420.0
        am(6, 3) = -AL * AL / 140.0
        am(4, 4) = 1.0 / 3.0
        am(5, 4) = 0.0
        am(6, 4) = 0.0
        am(5, 5) = 13.0 / 35.0
        am(6, 5) = -11.0 * AL / 210.0
        am(6, 6) = AL * AL / 105.0
        '対称項セット
        For i = 1 To 5
            For j = i + 1 To 6
                am(i, j) = am(j, i)
            Next j
        Next i
        For i = 1 To 6
            For j = 1 To 6
                am(i, j) = gamma(ne) * AA(ne) * AL * am(i, j) / ggg
            Next j
        Next i

    End Sub


    Private Sub GLTMX(ByVal ne As Integer, ByRef kakom(,) As Integer, ByRef x() As Double, ByRef y() As Double, _
                     ByRef gt(,) As Double, ByRef hen(,) As Double)
        '座標変換行列と全体座標系への行列変換
        Dim i As Integer
        Dim j As Integer
        Dim k As Integer
        Dim w(6, 6) As Double
        Dim AL As Double
        Dim cs As Double
        Dim sn As Double

        '行列初期化
        For i = 1 To 6
            For j = 1 To 6
                hen(i, j) = 0.0
        Next j, i

        '要素長さ
        i = kakom(ne, 1)
        j = kakom(ne, 2)
        AL = Math.Sqrt((x(j) - x(i)) ^ 2 + (y(j) - y(i)) ^ 2)

        '座標変換行列
        cs = (x(j) - x(i)) / AL
        sn = (y(j) - y(i)) / AL
        hen(1, 1) = cs
        hen(1, 2) = sn
        hen(2, 1) = -sn
        hen(2, 2) = cs
        hen(3, 3) = 1.0
        For i = 1 To 3
            For j = 1 To 3
                hen(i + 3, j + 3) = hen(i, j)
        Next j, i

        '全体座標系への変換
        '(w)=(T)'*(gt)
        For i = 1 To 6
            For j = 1 To 6
                w(i, j) = 0.0
                For k = 1 To 6
                    w(i, j) = w(i, j) + hen(k, i) * gt(k, j)
        Next k, j, i
        '(w)*(T)=(T)'(gt)(T)：(sm)は全体座標系での要素剛性行列
        For i = 1 To 6
            For j = 1 To 6
                gt(i, j) = 0.0
                For k = 1 To 6
                    gt(i, j) = gt(i, j) + w(i, k) * hen(k, j)
        Next k, j, i

    End Sub


    Private Sub MV_FRM(ByVal ne As Integer, ByRef am(,) As Double, ByVal ggg As Double, _
        ByRef gkh() As Double, ByRef gkv() As Double, ByRef fevec() As Double)
        '慣性力ベクトル
        Dim i As Integer
        Dim j As Integer
        Dim work(6) As Double

        work(1) = ggg * gkh(ne)
        work(2) = ggg * gkv(ne)
        work(3) = 0.0
        work(4) = ggg * gkh(ne)
        work(5) = ggg * gkv(ne)
        work(6) = 0.0

        '全体座標系での慣性力ベクトル
        For i = 1 To 6
            fevec(i) = 0.0
            For j = 1 To 6
                fevec(i) = fevec(i) + am(i, j) * work(j)
            Next j
        Next i
    End Sub


    Private Sub TM_FRM(ByVal ne As Integer, ByRef kakom(,) As Integer, ByRef deltaT() As Double, _
        ByRef Em() As Double, ByRef AA() As Double, ByRef alpha() As Double, ByRef fevec() As Double, ByRef hen(,) As Double)
        '温度荷重相当ベクトル
        Dim i As Integer
        Dim j As Integer
        Dim tempa As Double
        Dim work(6) As Double

        i = kakom(ne, 1)
        j = kakom(ne, 2)
        tempa = 0.5 * (deltaT(i) + deltaT(j))
        '要素座標系での温度荷重相当ベクトル
        work(1) = -1.0 * Em(ne) * AA(ne) * alpha(ne) * tempa
        work(2) = 0.0
        work(3) = 0.0
        work(4) = 1.0 * Em(ne) * AA(ne) * alpha(ne) * tempa
        work(5) = 0.0
        work(6) = 0.0
        '全体座標系での温度荷重相当ベクトル
        '（座標変換行列は転置=逆行列を要素座標系でのベクトルに乗じる）
        For i = 1 To 6
            fevec(i) = 0.0
            For j = 1 To 6
                fevec(i) = fevec(i) + hen(j, i) * work(j)
            Next j
        Next i
    End Sub


    Private Sub CALST(ByVal ne As Integer, ByRef kakom(,) As Integer, ByRef dis() As Double, _
                      ByRef sm(,) As Double, ByRef hen(,) As Double, ByRef reac() As Double, ByRef tsm(,) As Double, _
                      ByRef Em() As Double, ByRef AA() As Double, ByRef alpha() As Double, ByRef deltaT() As Double)
        '断面力算定
        Dim i As Integer
        Dim j As Integer
        Dim k As Integer
        Dim wd(6) As Double
        Dim ws(6) As Double
        Dim ia As Integer
        Dim ir As Integer
        Dim iw As Integer
        Dim tempa As Double

        '要素節点変位取得（全体座標系）
        For i = 1 To 2
            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)
        Next i
        '全体座標系での要素節点力算定
        For i = 1 To 6
            ws(i) = 0.0
            For k = 1 To 6
                ws(i) = ws(i) + sm(i, k) * wd(k)
            Next k
        Next i
        '要素断面力算定：座標変換行列×全体座標系要素節点力
        '今後使用しない全体剛性行列格納配列(tsm)に断面力を格納
        For i = 1 To 6
            tsm(ne, i) = 0.0
            For k = 1 To 6
                tsm(ne, i) = tsm(ne, i) + hen(i, k) * ws(k)
            Next k
        Next i
        '軸力の温度応力補正
        i = kakom(ne, 1)
        j = kakom(ne, 2)
        tempa = 0.5 * (deltaT(i) + deltaT(j))
        tsm(ne, 1) = tsm(ne, 1) + Em(ne) * AA(ne) * alpha(ne) * tempa
        tsm(ne, 4) = tsm(ne, 4) - Em(ne) * AA(ne) * alpha(ne) * tempa
        '全体座標系での内力ベクトル
        '（座標変換行列は転置=逆行列を要素座標系でのベクトルに乗じる）
        For i = 1 To 6
            ws(i) = 0.0
            For j = 1 To 6
                ws(i) = ws(i) + hen(j, i) * tsm(ne, j)
            Next j
        Next i

        '要素節点力の全体系へのセット
        For i = 1 To 2
            ia = kakom(ne, i)
            For j = 1 To 3
                ir = 3 * (ia - 1) + j
                iw = 3 * (i - 1) + j
                reac(ir) = reac(ir) + ws(iw)
            Next j
        Next i

    End Sub


    Private Sub ASMAT(ByVal ne As Integer, ByVal nod As Integer, ByVal nhen As Integer, _
                      ByRef kakom(,) As Integer, ByRef sm(,) As Double, ByRef tsm(,) As Double)
        '剛性行列組立(nod：1要素節点数，nhen：1節点の自由度)
        Dim i As Integer : Dim j As Integer : Dim k As Integer : Dim l As Integer
        Dim ki As Integer : Dim kj As Integer : Dim ks As Integer : Dim js As Integer
        Dim kik As Integer : Dim kjl As Integer : Dim ksk As Integer : Dim jsl As Integer

        For i = 1 To nod
            For j = 1 To nod
                ki = (kakom(ne, i) - 1) * nhen
                kj = (kakom(ne, j) - 1) * nhen
                ks = (i - 1) * nhen
                js = (j - 1) * nhen
                For k = 1 To nhen
                    For l = 1 To nhen
                        kik = ki + k : kjl = kj + l : ksk = ks + k : jsl = js + l
                        tsm(kik, kjl) = tsm(kik, kjl) + sm(ksk, jsl)
                Next l, k
        Next j, i

    End Sub


    Private Sub ASVEC(ByVal ne As Integer, ByVal nod As Integer, ByVal nhen As Integer, _
               ByRef kakom(,) As Integer, ByRef fevec() As Double, ByRef ftvec() As Double)
        '熱流束ベクトル組立(nod：1要素節点数，nhen：1節点の自由度)
        Dim i As Integer : Dim k As Integer
        Dim ki As Integer : Dim ks As Integer : Dim kik As Integer : Dim ksk As Integer

        For i = 1 To nod
            ki = (kakom(ne, i) - 1) * nhen
            ks = (i - 1) * nhen
            For k = 1 To nhen
                kik = ki + k : ksk = ks + k
                ftvec(kik) = ftvec(kik) + fevec(ksk)
            Next k
        Next i

    End Sub


    Private Sub CHBAND10(ByVal n1 As Integer, ByVal m1 As Integer, ByRef array(,) As Double)
        Dim i As Integer
        Dim j As Integer
        Dim k As Integer
        Dim mm As Integer
        Dim mn As Integer
        Dim z As Double

        array(1, 1) = Math.Sqrt(array(1, 1))
        For j = 2 To m1
            array(1, j) = array(1, j) / array(1, 1)
        Next j
        For i = 2 To n1
            If n1 - i + 1 > m1 Then
                mm = m1
            Else
                mm = n1 - i + 1
            End If
            For j = 1 To mm
                z = array(i, j)
                If j <> m1 Then
                    If m1 - j + 1 < i Then
                        mn = m1 - j + 1
                    Else
                        mn = i
                    End If
                    For k = 2 To mn
                        z = z - array(i - k + 1, k) * array(i - k + 1, j + k - 1)
                    Next k
                    If j > 1 Then
                        array(i, j) = z / array(i, 1)
                    Else
                        array(i, 1) = Math.Sqrt(z)
                    End If
                Else
                    array(i, j) = z / array(i, 1)
                End If
        Next j, i

    End Sub


    Private Sub CHBAND11(ByVal n1 As Integer, ByVal m1 As Integer, ByRef array(,) As Double, ByRef vector() As Double)
        Dim i As Integer
        Dim j As Integer
        Dim mm As Integer
        Dim z As Double

        vector(1) = vector(1) / array(1, 1)
        For i = 2 To n1
            z = vector(i)
            If i > m1 Then
                mm = m1
            Else
                mm = i
            End If
            For j = 2 To mm
                z = z - array(i - j + 1, j) * vector(i - j + 1)
            Next j
            vector(i) = z / array(i, 1)
        Next i

        vector(n1) = vector(n1) / array(n1, 1)
        For i = n1 - 1 To 1 Step -1
            z = vector(i)
            For j = 2 To m1
                If i + j - 1 > n1 Then Exit For
                z = z - array(i, j) * vector(i + j - 1)
            Next j
            vector(i) = z / array(i, 1)
        Next i

    End Sub


    Private Function IBAND(ByVal mm As Integer, ByRef array(,) As Double, ByVal fact0 As Double) As Integer
        Dim i As Integer
        Dim j As Integer
        '-----------------------------------------------------
        'バンド幅計算
        '-----------------------------------------------------
        IBAND = 1
        For i = 1 To mm - 1
            For j = mm To i + 1 Step -1
                If fact0 <= Math.Abs(array(i, j)) Then Exit For
            Next j
            If IBAND <= j - i + 1 Then IBAND = j - i + 1
        Next i
        If mm < IBAND Then IBAND = mm

    End Function


    Private Sub BAND(ByVal mm As Integer, ByVal ib As Integer, ByRef array(,) As Double)
        Dim i As Integer
        Dim j As Integer
        Dim mn As Integer
        Dim work() As Double
        ReDim work(mm)
        '-----------------------------------------------------
        'フルマトリックスのバンド化記憶
        '-----------------------------------------------------
        For i = 2 To mm
            If i + ib - 1 < mm Then
                mn = ib
            Else
                mn = mm - i + 1
            End If
            For j = 1 To ib
                work(j) = 0.0
            Next j
            For k = 1 To mn
                work(k) = array(i, i + k - 1)
            Next k
            For j = 1 To ib
                array(i, j) = work(j)
            Next j
        Next i

    End Sub


    Private Sub IERROR2(ByVal fnameW As String)
        '結果出力（Appendモード）
        Dim sw As System.IO.StreamWriter
        Dim dat As String
        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", "異常終了通知")
        Me.Close()

    End Sub

End Class