constructcls.cls

该程序是按照矩阵位移法的后处理法的基本原理和分析过程

    ' WRITE NODAL DISPLACEMENTS
    '写结点位移
    Print #fp, "NODAL DISPLACEMENTS"
    Print #fp, Tab(8); "NODE"; Tab(30); "U"; Tab(50); "V"
    For i = 1 To NN
        K1 = NDF * (i - 1) + 1
        K2 = K1 + NDF - 1
        Print #fp, Tab(10); Format(i, "00");
        For j = K1 To K2
           Print #fp, Tab(10 + 20 * (j - K1 + 1)); Format(AL(j), "0.000000");
        Next j
        Print #fp, " "
    Next i
    Print #fp, ""
' WRITE NODAL REACTIONS
    Print #fp, "NODAL REACTIONS"
    Print #fp, Tab(10); "NODE"; Tab(30); "PX"; Tab(50); "PY"
    For i = 1 To NN
        K1 = NDF * (i - 1) + 1
        K2 = K1 + NDF - 1
        Print #fp, Spc(10); Format(i, "00");
        For j = K1 To K2
            Print #fp, Tab(10 + 20 * (j - K1 + 1)); Format(REAC(j), "0.000000");
        Next j
        Print #fp, " "
    Next i
    Print #fp, ""
'WRITE MEMBER AXIAL FORCES
    Print #fp, "MEMBER FORCES"; "      MEMBER    AXIAL FORCE"
    For i = 1 To NE
       Print #fp, Tab(10); Format(i, "00"); Tab(20); Format(FORC(i), "0.000000")
    Next i
    Print #fp, ""
    Print #fp, ""
    Print #fp, "**********************************************************************";
    Close #fp
End Sub
'****************************************************************************************

'结果形成和计算过程
'*****************************************************************************************
'计算
Public Sub calculate()
    
    If HaveReaded = False Then Exit Sub
    ModifyIB
    ASSEM
    BOUND
    SLBSI
    Force
    GetXYDis 1, scale1
    isCalculate = True
    ModifyIBback
End Sub

Private Sub ModifyIB()
    Dim h As Integer
    For h = 1 To NBN
        If IB(3 * h - 1) <> 0 Then
            IB(3 * h - 1) = 0
        Else: IB(3 * h - 1) = 1
        End If
        If IB(3 * h) <> 0 Then
            IB(3 * h) = 0
        Else: IB(3 * h) = 1
        End If
    Next h
End Sub
'
Private Sub ModifyIBback()
    IB = IB1
End Sub

'形成总刚度,确定矩阵尺度界限NEMX,NCMX
Private Sub ASSEM()
    Dim N1 As Integer, i As Integer, L1 As Integer, j As Integer, _
        L2 As Integer, J1 As Integer, K As Integer, L3 As Integer, L As Integer, NEL As Integer
'COMPUTE HALF BAND WIDTH AND STORE IN MS
    N1 = NNE - 1
    MS = 0
    For i = 1 To NE
        L1 = NNE * (i - 1)
        For j = 1 To N1
            L2 = L1 + j
            J1 = j + 1
            For K = J1 To NNE
                L3 = L1 + K
                L = Abs(NCO(L2) - NCO(L3))
                If (MS - L) <= 0 Then MS = L
            Next K
        Next j
    Next i
    MS = NDF * (MS + 1)

'CLEAR THE TOTAL STIFFNESS MATRIX
    ReDim TK(N, MS) As Double
    ReDim ELST(4, 4) As Double
    For NEL = 1 To NE
'      COMPUTE THE STIFFNESS MATRIX FOR ELEMENT NEL
        '形成单刚
        Call STIFF(NEL)
'      PLACE THE MATRIX IN THE TOTAL STIFFNESS MATRIX
        '单刚送总刚
        Call ELASS(NEL, TK, ELST)
    Next NEL
End Sub

'形成单刚
' COMPUTATION OF ELEMENT STIFFNESS MATRIX FOR CURRENT ELEMENT
Private Sub STIFF(ByVal NEL As Integer)
    Dim L As Integer, N1 As Integer, N2 As Integer, i As Integer, j As Integer, _
        K1 As Integer, K2 As Integer
    Dim D As Double, CO As Double, SI As Double, COEF As Double
    L = NNE * (NEL - 1)
    N1 = NCO(L + 1)
    N2 = NCO(L + 2)
' COMPUTE LENGTH OF ELEMENT, AND SINE AND COSINE OF ITS LOCAL X AXIS
    D = Sqr((X(N2) - X(N1)) ^ 2 + (Y(N2) - Y(N1)) ^ 2)
    CO = (X(N2) - X(N1)) / D
    SI = (Y(N2) - Y(N1)) / D
' COMPUTE ELEMENT STIFFNESS MATRIX
    COEF = E * PROP(NEL) / D
    ELST(1, 1) = COEF * CO * CO
    ELST(1, 2) = COEF * CO * SI
    ELST(2, 2) = COEF * SI * SI
    For i = 1 To 2
       For j = i To 2
           K1 = i + NDF
           K2 = j + NDF
           ELST(K1, K2) = ELST(i, j)
           ELST(i, K2) = -ELST(i, j)
        Next j
    Next i
    ELST(2, 3) = -COEF * SI * CO
End Sub

'单刚送总刚
'STORE THE ELEMENT MATRIX FOR ELEMENT NEL IN THE TOTAL MATRIX
Private Sub ELASS(NEL As Integer, TM() As Double, ELMAT() As Double)
    Dim L1%, i%, L2%, N1%, I1%, J1%, j%, N2%, I2%, J2%, K%, KI%, KR%, IC%, K1%, K2%, L%, KC%
    L1 = NNE * (NEL - 1)
    For i = 1 To NNE
        L2 = L1 + i
        N1 = NCO(L2)
        I1 = NDF * (i - 1)
        J1 = NDF * (N1 - 1)
        For j = i To NNE
            L2 = L1 + j
            N2 = NCO(L2)
            I2 = NDF * (j - 1)
            J2 = NDF * (N2 - 1)
            For K = 1 To NDF
                KI = 1
                If (N1 - N2) = 0 Then
                    KI = K        'STORE A DIAGONAL SUBMATRIX
                End If
                If (N1 - N2) <= 0 Then              ' STORE AN OFF DIAGONAL SUBMATRIX
                    KR = J1 + K
                    IC = J2 - KR + 1
                    K1 = I1 + K
                Else                'STORE THE TRANSPOSE OF AN OFF DIAGONAL MATRIX
                    KR = J2 + K
                    IC = J1 - KR + 1
                    K2 = I2 + K
                End If
                For L = KI To NDF
                    KC = IC + L
                    If (N1 - N2) <= 0 Then
                        K2 = I2 + L
                    Else
                        K1 = I1 + L
                    End If
                    TM(KR, KC) = TM(KR, KC) + ELMAT(K1, K2)
                Next L
            Next K
        Next j
    Next i
End Sub

'边界条件处理
'INTRODUCTION OF THE BOUNDARY CONDITIONS
Private Sub BOUND()
    Dim L%, L1%, NO%, K1%, i%, L2%, KR%, j%, KV%
    For L = 1 To NBN
        L1 = (NDF + 1) * (L - 1) + 1
        NO = IB(L1)
        K1 = NDF * (NO - 1)
        For i = 1 To NDF
            L2 = L1 + i
            If IB(L2) = 0 Then
'              PRESCRIBED UNKNOWN TO BE CONSIDERED
                KR = K1 + i
                For j = 2 To MS
                    KV = KR + j - 1
                    If (N - KV) >= 0 Then
'                      MODIFY ROW OF TK AND CORRESPONDINF ELEMENTS IN AL
                        AL(KV) = AL(KV) - TK(KR, j) * REAC(KR)
                        TK(KR, j) = 0#
                    End If
                    KV = KR - j + 1
                    If KV > 0 Then
'                      MODIFY COLUMN IN TK AND CORRESPONDING ELEMENTS IN AL
                        AL(KV) = AL(KV) - TK(KV, j) * REAC(KR)
                        TK(KV, j) = 0#
                    End If
                Next j
'              SET DIAGONAL COEFFICIENT OF TK EQUAL TO 1 PLACE PRESCRIBED UNKNOWN
'              VALUE IN AL
                TK(KR, 1) = 1#
                AL(KR) = REAC(KR)
            End If
        Next i
    Next L
End Sub

' SOLUTION OF SIMUTANEOUS SYSTEMS OF EQUATIONS BY THE GAUSS
' ELIMINATION METHOD,FOR SYMMETRIC BANDED MATRICES
'求解方程组        TK           AL              V
Private Sub SLBSI()
    Dim N1%, K%, K1%, NI%, L%, j%, K2%, i%, K3%
    Dim C As Double
    ReDim V(20) As Double
    N1 = N - 1
    For K = 1 To N1
        C = TK(K, 1)
        K1 = K + 1
        Dim s As String
        If C <= 0.000001 And C >= -0.000001 Then
            s = "  **** SINGULARITY IN ROW " + Str$(K)
            MsgBox s
            Exit Sub
        Else                                'DIVIDE ROW BY DIAGONAL COEFFICIENT
            NI = K1 + MS - 2
            If NI <= N Then L = NI
            If NI > N Then L = N
            For j = 2 To MS
                V(j) = TK(K, j)
            Next j
            For j = K1 To L
                K2 = j - K + 1
                TK(K, K2) = TK(K, K2) / C
            Next j
            AL(K) = AL(K) / C
'          ELIMINATE UNKNOWN X(K) FROM ROW I
            For i = K1 To L
                K2 = i - K1 + 2
                C = V(K2)
                For j = i To L
                    K2 = j - i + 1
                    K3 = j - K + 1
                    TK(i, K2) = TK(i, K2) - C * TK(K, K3)
                Next j
                AL(i) = AL(i) - C * AL(K)
            Next i
        End If
    Next K
'   COMPUTE LAST UNKNOWN
    If TK(N, 1) <= 0.000001 And TK(N, 1) >= 0.000001 Then
            s = "  **** SINGULARITY IN ROW " + Str$(K)
            MsgBox s
        Exit Sub
    Else
        AL(N) = AL(N) / TK(N, 1)
'      APPLY BACKSUBSTITUTE PROCESS TO COMPUTE REMAINING UNKNOWNS
        For i = 1 To N1
            K = N - i
            K1 = K + 1
            NI = K1 + MS - 2
            If NI <= N Then L = NI
            If NI > N Then L = N
            For j = K1 To L
                K2 = j - K + 1
                AL(K) = AL(K) - TK(K, K2) * AL(j)
            Next j
        Next i
    End If
End Sub

' COMPUTATION OF ELEMENT FORCES
'计算杆端力
Private Sub Force()
    Dim i%, NEL%, L%, N1%, N2%, K1%, K2%
    Dim D As Double, CO As Double, SI As Double, COEF As Double
' CLEAR THE REACTIONS ARRAY
    ReDim REAC(N) As Double
    ReDim FORC(NE) As Double
    For NEL = 1 To NE
        L = NNE * (NEL - 1)
        N1 = NCO(L + 1)
        N2 = NCO(L + 2)
        K1 = NDF * (N1 - 1)
        K2 = NDF * (N2 - 1)
'      COMPUTE LENGTH OF ELEMENT, AND SINE/COSINE OF ITS LOCAL X AXIS
        D = Sqr((X(N2) - X(N1)) ^ 2 + (Y(N2) - Y(N1)) ^ 2)
        CO = (X(N2) - X(N1)) / D
        SI = (Y(N2) - Y(N1)) / D
        COEF = E * PROP(NEL) / D
'      COMPUTE MEMBER AXIAL FORCE AND STORE IN ARRAY FORC
        FORC(NEL) = COEF * ((AL(K2 + 1) - AL(K1 + 1)) * CO + (AL(K2 + 2) - AL(K1 + 2)) * SI)
'      COMPUTE NODAL RESULTANTS
        REAC(K1 + 1) = REAC(K1 + 1) - FORC(NEL) * CO
        REAC(K1 + 2) = REAC(K1 + 2) - FORC(NEL) * SI
        REAC(K2 + 1) = REAC(K2 + 1) + FORC(NEL) * CO
        REAC(K2 + 2) = REAC(K2 + 2) + FORC(NEL) * SI
    Next NEL
End Sub
'**********************************************************************************************
Private Sub Class_Terminate()

End Sub