Option Explicit On
Option Strict On
Public Class Form1

    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 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 fname1 As String = ""
        Dim fname2 As String = ""
        Dim k As Integer
        Dim iii As Integer
        Dim jjj As Integer
        Dim strcom As String 'コメント文
        Dim nd As Integer '計算組数
        Dim STRES As Integer '0:平面歪み，1：平面応力
        Dim P() As Double '内圧
        Dim a() As Double '水路内半径
        Dim b() As Double '掘削半径
        Dim Ec() As Double 'コンクリート弾性係数
        Dim nc() As Double 'コンクリートポアソン比
        Dim Es() As Double '鉄筋弾性係数
        Dim ta() As Double '内側鉄筋等価板厚
        Dim tb() As Double '外側鉄筋等価板厚
        Dim Eg() As Double '岩盤弾性係数
        Dim ng() As Double '岩盤ポアソン比
        Dim sigsa(,) As Double '内側鉄筋応力
        Dim sigsb(,) As Double '外側鉄筋応力
        Dim sigrc(,) As Double 'コンクリート半径方向応力
        Dim sigtc(,) As Double 'コンクリート円周方向応力
        Dim ua(,) As Double '内壁半径方向変位
        Dim ub(,) As Double '外壁半径方向変位
        '作業用変数
        Dim PP As Double : Dim ar As Double : Dim br As Double
        Dim Ecc As Double : Dim ncc As Double
        Dim Ess As Double : Dim taa As Double : Dim tbb As Double
        Dim Egg As Double : Dim ngg As Double
        Dim ssa As Double : Dim ssb As Double
        Dim src As Double : Dim stc As Double
        Dim uaa As Double : Dim ubb As Double

        '入力ファイル選定とテキスト入力
        If OpenFileDialog1.ShowDialog() = Windows.Forms.DialogResult.OK Then fname1 = OpenFileDialog1.FileName
        sr = New System.IO.StreamReader(fname1, System.Text.Encoding.Default)
        dat = sr.ReadLine() : sbuf = dat.Split(delim) : strcom = sbuf(0)
        dat = sr.ReadLine() : sbuf = dat.Split(delim) : nd = CInt(sbuf(0)) - 1 : STRES = CInt(sbuf(1))
        ReDim P(nd) : ReDim a(nd) : ReDim b(nd)
        ReDim Ec(nd) : ReDim nc(nd)
        ReDim Es(nd) : ReDim ta(nd) : ReDim tb(nd)
        ReDim Eg(nd) : ReDim ng(nd)
        k = 0
        Do While sr.Peek() >= 0
            dat = sr.ReadLine() : sbuf = dat.Split(delim)
            P(k) = CDbl(sbuf(0))
            a(k) = CDbl(sbuf(1))
            b(k) = CDbl(sbuf(2))
            Ec(k) = CDbl(sbuf(3))
            nc(k) = CDbl(sbuf(4))
            Es(k) = CDbl(sbuf(5))
            ta(k) = CDbl(sbuf(6))
            tb(k) = CDbl(sbuf(7))
            Eg(k) = CDbl(sbuf(8))
            ng(k) = CDbl(sbuf(9))
            k = k + 1
        Loop
        sr.Close()
        nd = k - 1

        ReDim sigsa(nd, 3)
        ReDim sigsb(nd, 3)
        ReDim sigrc(nd, 3)
        ReDim sigtc(nd, 3)
        ReDim ua(nd, 3)
        ReDim ub(nd, 3)
        For iii = 0 To nd
            PP = P(iii) : ar = a(iii) : br = b(iii)
            Ecc = Ec(iii) : ncc = nc(iii)
            Ess = Es(iii) : taa = ta(iii) : tbb = tb(iii)
            Egg = Eg(iii) : ngg = ng(iii)
            For jjj = 0 To 3
                Select Case jjj
                    Case 0 'ひび割れのない無筋コンクリート
                        Call CALC0(STRES, PP, ar, br, Ecc, ncc, Ess, taa, tbb, Egg, ngg, ssa, ssb, src, stc, uaa, ubb)
                    Case 1 'ひび割れの無い単鉄筋コンクリート
                        Call CALC1(STRES, PP, ar, br, Ecc, ncc, Ess, taa, tbb, Egg, ngg, ssa, ssb, src, stc, uaa, ubb)
                    Case 2 'ひび割れのある単鉄筋コンクリート
                        Call CALC2(STRES, PP, ar, br, Ecc, ncc, Ess, taa, tbb, Egg, ngg, ssa, ssb, src, stc, uaa, ubb)
                    Case 3 'ひび割れのある複鉄筋コンクリート
                        Call CALC3(STRES, PP, ar, br, Ecc, ncc, Ess, taa, tbb, Egg, ngg, ssa, ssb, src, stc, uaa, ubb)
                End Select
                sigsa(iii, jjj) = ssa '内側鉄筋応力
                sigsb(iii, jjj) = ssb '外側鉄筋応力
                sigrc(iii, jjj) = src 'コンクリート半径方向応力
                sigtc(iii, jjj) = stc 'コンクリート円周方向応力
                ua(iii, jjj) = uaa '内壁半径方向変位
                ub(iii, jjj) = ubb '外壁半径方向変位
            Next jjj
        Next iii

        '出力ファイル選定とテキスト出力
        If SaveFileDialog1.ShowDialog() = Windows.Forms.DialogResult.OK Then fname2 = SaveFileDialog1.FileName
        sw = New System.IO.StreamWriter(fname2, False, System.Text.Encoding.Default)
        dat = strcom : sw.WriteLine(dat)
        dat = "nd," & (nd + 1).ToString & ",STRES," & STRES.ToString : sw.WriteLine(dat)
        dat = "入力データ" : sw.WriteLine(dat)
        dat = "i,P,a,b,Ec,nc,Es,ta,tb,Eg,ng" : sw.WriteLine(dat)
        For iii = 0 To nd
            dat = CStr(iii + 1) & "," & P(iii) & "," & a(iii) & "," & b(iii) & "," & Ec(iii) & "," & nc(iii) & "," & Es(iii) & "," & ta(iii) & "," & tb(iii) & "," & Eg(iii) & "," & ng(iii)
            sw.WriteLine(dat)
        Next iii
        For jjj = 0 To 3
            Select Case jjj
                Case 0 : dat = "ひび割れのない無筋コンクリート"
                Case 1 : dat = "ひび割れのない単鉄筋コンクリート"
                Case 2 : dat = "ひび割れのある単鉄筋コンクリート"
                Case 3 : dat = "ひび割れのある複鉄筋コンクリート"
            End Select
            sw.WriteLine(dat)
            dat = "i,σsa,σsb,σrc,σtc,u(r=a),u(r=b)" : sw.WriteLine(dat)
            For iii = 0 To nd
                dat = CStr(iii + 1) & "," & sigsa(iii, jjj) & "," & sigsb(iii, jjj) & "," & sigrc(iii, jjj) & "," & sigtc(iii, jjj) & "," & ua(iii, jjj) & "," & ub(iii, jjj)
                sw.WriteLine(dat)
            Next iii
        Next jjj
        sw.Close()

    End Sub

    Private Sub STIF(ByVal STRES As Integer, ByVal EE As Double, ByVal po As Double, ByRef AA As Double, ByRef BB As Double)
        Select Case STRES
            Case 0 '平面ひずみ
                AA = EE * (1.0 - po) / (1.0 + po) / (1.0 - 2.0 * po)
                BB = EE * po / (1.0 + po) / (1.0 - 2.0 * po)
            Case 1 '平面応力
                AA = EE / (1.0 - po * po)
                BB = EE * po / (1.0 - po * po)
        End Select
    End Sub

    Private Sub CALC0(ByRef STRES As Integer, ByVal PP As Double, ByVal ar As Double, ByVal br As Double, _
                      ByVal Ecc As Double, ByVal ncc As Double, _
                      ByVal Ess As Double, ByVal taa As Double, ByVal tbb As Double, _
                      ByVal Egg As Double, ByVal ngg As Double, _
                      ByRef ssa As Double, ByRef ssb As Double, _
                      ByRef src As Double, ByRef stc As Double, _
                      ByRef uaa As Double, ByRef ubb As Double)
        'ひび割れのない無筋コンクリート
        Dim sm(,) As Double : Dim mm As Integer
        Dim Ac As Double : Dim Bc As Double
        Dim C1 As Double : Dim C2 As Double : Dim Cg As Double

        mm = 2
        Call STIF(STRES, Ecc, ncc, Ac, Bc)
        ReDim sm(mm, mm + 1)
        For i = 0 To mm
            For j = 0 To mm + 1
                sm(i, j) = 0.0
        Next j, i
        sm(0, 0) = Ac + Bc
        sm(0, 1) = -(Ac - Bc) / ar ^ 2
        sm(1, 0) = br
        sm(1, 1) = 1.0 / br
        sm(1, 2) = -1.0 / br
        sm(2, 0) = Ac + Bc
        sm(2, 1) = -(Ac - Bc) / br ^ 2
        sm(2, 2) = Egg / (1.0 + ngg) / br ^ 2
        sm(0, 3) = -PP
        Call MATGJ(mm, sm)
        C1 = sm(0, 3)
        C2 = sm(1, 3)
        Cg = sm(2, 3)
        ssa = 0.0
        ssb = 0.0
        src = (Ac + Bc) * C1 - (Ac - Bc) * C2 / ar ^ 2
        stc = (Ac + Bc) * C1 + (Ac - Bc) * C2 / ar ^ 2
        uaa = C1 * ar + C2 / ar
        ubb = C1 * br + C2 / br

    End Sub

    Private Sub CALC1(ByRef STRES As Integer, ByVal PP As Double, ByVal ar As Double, ByVal br As Double, _
                       ByVal Ecc As Double, ByVal ncc As Double, _
                       ByVal Ess As Double, ByVal taa As Double, ByVal tbb As Double, _
                       ByVal Egg As Double, ByVal ngg As Double, _
                       ByRef ssa As Double, ByRef ssb As Double, _
                       ByRef src As Double, ByRef stc As Double, _
                       ByRef uaa As Double, ByRef ubb As Double)
        'ひび割れのない単鉄筋コンクリート
        Dim sm(,) As Double : Dim mm As Integer
        Dim Ac As Double : Dim Bc As Double
        Dim C1 As Double : Dim C2 As Double : Dim Cg As Double : Dim Psa As Double

        mm = 3
        Call STIF(STRES, Ecc, ncc, Ac, Bc)
        ReDim sm(mm, mm + 1)
        For i = 0 To mm
            For j = 0 To mm + 1
                sm(i, j) = 0.0
        Next j, i
        sm(0, 0) = ar
        sm(0, 1) = 1.0 / ar
        sm(0, 2) = -ar ^ 2 / Ess / taa
        sm(1, 0) = Ac + Bc
        sm(1, 1) = -(Ac - Bc) / ar ^ 2
        sm(1, 2) = -1.0
        sm(2, 0) = br
        sm(2, 1) = 1.0 / br
        sm(2, 3) = -1.0 / br
        sm(3, 0) = Ac + Bc
        sm(3, 1) = -(Ac - Bc) / br ^ 2
        sm(3, 3) = Egg / (1.0 + ngg) / br ^ 2
        sm(1, 4) = -PP
        Call MATGJ(mm, sm)
        C1 = sm(0, 4)
        C2 = sm(1, 4)
        Psa = sm(2, 4)
        Cg = sm(3, 4)
        ssa = Psa * ar / taa
        ssb = 0.0
        src = (Ac + Bc) * C1 - (Ac - Bc) * C2 / ar ^ 2
        stc = (Ac + Bc) * C1 + (Ac - Bc) * C2 / ar ^ 2
        uaa = Psa * ar ^ 2 / Ess / taa
        ubb = C1 * br + C2 / br

    End Sub

    Private Sub CALC2(ByRef STRES As Integer, ByVal PP As Double, ByVal ar As Double, ByVal br As Double, _
                    ByVal Ecc As Double, ByVal ncc As Double, _
                    ByVal Ess As Double, ByVal taa As Double, ByVal tbb As Double, _
                    ByVal Egg As Double, ByVal ngg As Double, _
                    ByRef ssa As Double, ByRef ssb As Double, _
                    ByRef src As Double, ByRef stc As Double, _
                    ByRef uaa As Double, ByRef ubb As Double)
        'ひび割れのある単鉄筋コンクリート
        Dim sm(,) As Double : Dim mm As Integer
        Dim C3 As Double : Dim C4 As Double : Dim Cg As Double : Dim Psa As Double

        mm = 3
        ReDim sm(mm, mm + 1)
        For i = 0 To mm
            For j = 0 To mm + 1
                sm(i, j) = 0.0
        Next j, i
        sm(0, 0) = Math.Log(ar) / Ecc
        sm(0, 1) = 1.0
        sm(0, 2) = -ar ^ 2 / Ess / taa
        sm(1, 0) = 1.0 / ar
        sm(1, 2) = -1.0
        sm(2, 0) = Math.Log(br) / Ecc
        sm(2, 1) = 1.0
        sm(2, 3) = -1.0 / br
        sm(3, 0) = 1.0 / br
        sm(3, 3) = Egg / (1.0 + ngg) / br ^ 2
        sm(1, 4) = -PP
        Call MATGJ(mm, sm)
        C3 = sm(0, 4)
        C4 = sm(1, 4)
        Psa = sm(2, 4)
        Cg = sm(3, 4)
        ssa = Psa * ar / taa
        ssb = 0.0
        src = C3 / ar
        stc = 0.0
        uaa = Psa * ar ^ 2 / Ess / taa
        ubb = C3 / Ecc * Math.Log(br) + C4

    End Sub

    Private Sub CALC3(ByRef STRES As Integer, ByVal PP As Double, ByVal ar As Double, ByVal br As Double, _
                     ByVal Ecc As Double, ByVal ncc As Double, _
                     ByVal Ess As Double, ByVal taa As Double, ByVal tbb As Double, _
                     ByVal Egg As Double, ByVal ngg As Double, _
                     ByRef ssa As Double, ByRef ssb As Double, _
                     ByRef src As Double, ByRef stc As Double, _
                     ByRef uaa As Double, ByRef ubb As Double)
        'ひび割れのある複鉄筋コンクリート
        Dim sm(,) As Double : Dim mm As Integer
        Dim C3 As Double : Dim C4 As Double : Dim Cg As Double : Dim Psa As Double : Dim Psb As Double

        If 0.0 < tbb Then
            mm = 4
            ReDim sm(mm, mm + 1)
            For i = 0 To mm
                For j = 0 To mm + 1
                    sm(i, j) = 0.0
            Next j, i
            sm(0, 0) = Math.Log(ar) / Ecc
            sm(0, 1) = 1.0
            sm(0, 2) = -ar ^ 2 / Ess / taa
            sm(1, 0) = 1.0 / ar
            sm(1, 2) = -1.0
            sm(2, 0) = Math.Log(br) / Ecc
            sm(2, 1) = 1.0
            sm(2, 3) = -br ^ 2 / Ess / tbb
            sm(3, 3) = br ^ 2 / Ess / tbb
            sm(3, 4) = -1.0 / br
            sm(4, 0) = 1.0 / br
            sm(4, 3) = 1.0
            sm(4, 4) = Egg / (1.0 + ngg) / br ^ 2
            sm(1, 5) = -PP
            Call MATGJ(mm, sm)
            C3 = sm(0, 5)
            C4 = sm(1, 5)
            Psa = sm(2, 5)
            Psb = sm(3, 5)
            Cg = sm(4, 5)
            ssa = Psa * ar / taa
            ssb = Psb * br / tbb
            src = C3 / ar
            stc = 0.0
            uaa = Psa * ar ^ 2 / Ess / taa
            ubb = C3 / Ecc * Math.Log(br) + C4
        Else
            ssa = 0.0
            ssb = 0.0
            src = 0.0
            stc = 0.0
            uaa = 0.0
            ubb = 0.0
        End If

    End Sub

    Private Sub MATGJ(ByVal n As Integer, ByRef a(,) As Double)
        'Gauss-Jordan法による連立一次方程式の解法
        Dim i As Integer
        Dim j As Integer
        Dim k As Integer
        Dim s As Integer
        Dim p As Double
        Dim d As Double
        Dim max As Double
        Dim dumy As Double

        For k = 0 To n
            max = 0.0
            s = k
            For j = k To n
                If Math.Abs(a(j, k)) > max Then
                    max = Math.Abs(a(j, k))
                    s = j
                End If
            Next j
            For j = 0 To n + 1
                dumy = a(k, j)
                a(k, j) = a(s, j)
                a(s, j) = dumy
            Next j
            p = a(k, k)
            For j = k To n + 1
                a(k, j) = a(k, j) / p
            Next j
            For i = 0 To n
                If i <> k Then
                    d = a(i, k)
                    For j = k To n + 1
                        a(i, j) = a(i, j) - d * a(k, j)
                    Next j
                End If
            Next i
        Next k

    End Sub

End Class