mdlfinite.bas

用visual basic编写的有限元程序！

Attribute VB_Name = "mdlFinite"
Private Const Max_Table = 1000

Public Const CN_NUL = -1
Public Const CN_X = 0
Public Const CN_Y = 1
Public Const CN_Z = 2
Public Const CN_XY = 3
Public Const CN_YZ = 4
Public Const CN_XZ = 5
Public Const CN_ALL = 6

Public Type Ans_Node
    Point(2) As Double
    Theta(2) As Double
    U_Code As Integer
    U(5) As Double   '移动,转动自由度
    F_Code As Integer
    F(5) As Double   '力,力矩
    Used As Long
    SEL As Integer
End Type

Public Type Ans_Node2
    U(5) As Double
    R(5) As Double   '约束反力
End Type

Public Type Ans_Elem
    MAT As Long
    TYP As Long
    REL As Long
    ESY As Long
    SEC As Long
    SEL As Integer
    Nodes(1 To 8) As Long
End Type

Public Type Ans_Elem2
    E As Double  '弹性模量
    U As Double  '泊松比
    t As Double  '厚度
    Area As Double  '面积
    a(1 To 3) As Double
    b(1 To 3) As Double
    C(1 To 3) As Double
    ArrD(1 To 3, 1 To 3) As Double
    ArrB(1 To 3, 1 To 6) As Double
    ArrK(1 To 6, 1 To 6) As Double
End Type

Public Type Ans_Mati
    Mode As Integer '0-各向同性;1-正交各向异性;2-各向异性
    Mode2 As Integer '0-Nu;1-Pr
    EX As Double
    EY As Double
    EZ As Double
    NUXY As Double
    NUYZ As Double
    NUXZ As Double
    GXY As Double
    GYZ As Double
    GXZ As Double
    ArrD(1 To 6, 1 To 6) As Double
End Type

Public Type Ans_Real
    R(1 To 6) As Double
End Type


Public DEF_ELEM As Ans_Elem  '单元默认设置
Public Nodes() As Ans_Node
Public Elems() As Ans_Elem

Public Matis() As Ans_Mati
Public Reals() As Ans_Real

Public Top_Node As Long
Public Top_Elem As Long

Public Top_Mati As Long
Public Top_Real As Long
Public Node2s() As Ans_Node2
Public Elem2s() As Ans_Elem2

Public Function ClearAll() As Long
    Rem 1.*** 清除内存
    ReDim Nodes(Max_Table), Node2s(Max_Table)
    ReDim Elems(Max_Table), Elem2s(Max_Table)
    ReDim Matis(Max_Table), Reals(Max_Table)
    Top_Node = 0: Top_Elem = 0
    Top_Mati = 0: Top_Real = 0
    With DEF_ELEM
        .TYP = 1
        .MAT = 1
        .REL = 1
        .ESY = 0
        .SEC = 1
    End With
End Function


Public Function Add_Node(Point#(), ByVal num&)
    Rem 添加节点
    If num = 0 Then
        num = Top_Node + 1
    End If
    Nodes(num).Point(0) = Point(0)
    Nodes(num).Point(1) = Point(1)
    If Top_Node < num Then Top_Node = num
End Function

Public Function Add_Elem(ByVal I&, ByVal J&, ByVal K&)
    Rem 添加单元
    Dim num&, Area#
    If num = 0 Then
        num = Top_Elem + 1
    End If
    Elems(num) = DEF_ELEM  '默认设置
    Area = Nodes(I).Point(0) * Nodes(J).Point(1) - Nodes(I).Point(1) * Nodes(J).Point(0)
    Area = Nodes(J).Point(0) * Nodes(K).Point(1) - Nodes(J).Point(1) * Nodes(K).Point(0) + Area
    Area = Nodes(K).Point(0) * Nodes(I).Point(1) - Nodes(K).Point(1) * Nodes(I).Point(0) + Area
    If Area > 0 Then
        Elems(num).Nodes(1) = I
        Elems(num).Nodes(2) = J
        Elems(num).Nodes(3) = K
    Else
        Elems(num).Nodes(1) = I
        Elems(num).Nodes(2) = K
        Elems(num).Nodes(3) = J
    End If
    If Top_Elem < num Then Top_Elem = num
End Function

Public Function Add_Dof(num&, sLab$, uValue#)
    Dim nCode As Integer
    Select Case UCase(sLab)
        Case "UX": nCode = CN_X
        Case "UY": nCode = CN_Y
        Case "UZ": nCode = CN_Z
        Case "ROTZ": nCode = CN_XY
        Case "ROTX": nCode = CN_YZ
        Case "ROTY": nCode = CN_XZ
        Case "ALL": nCode = CN_ALL
        Case Else: nCode = CN_NUL
    End Select
    If nCode >= CN_X And nCode <= CN_XZ Then
        Nodes(num).U_Code = Nodes(num).U_Code Or (2 ^ nCode)
        Nodes(num).U(nCode) = uValue
    ElseIf nCode = CN_ALL Then
        Nodes(num).U_Code = Nodes(num).U_Code Or 3
        Nodes(num).U(CN_X) = 0
        Nodes(num).U(CN_Y) = 0
    End If
End Function

Public Function Add_Force(num&, sLab$, uValue#)
    Dim nCode As Integer
    Select Case UCase(sLab)
        Case "FX": nCode = CN_X
        Case "FY": nCode = CN_Y
        Case "FZ": nCode = CN_Z
        Case "MZ": nCode = CN_XY
        Case "MX": nCode = CN_YZ
        Case "MY": nCode = CN_XZ
        Case "ALL": nCode = CN_ALL
        Case Else: nCode = CN_NUL
    End Select
    If nCode >= 0 Then
        Nodes(num).F_Code = Nodes(num).F_Code Or (2 ^ nCode)
        Nodes(num).F(nCode) = uValue
    End If
End Function
Public Function Add_MpData(num&, sLab$, uValue#)
    Dim nCode As Integer
    Select Case UCase(sLab)
        Case "EX": Matis(num).EX = uValue
        Case "EY": Matis(num).EY = uValue
        Case "EZ": Matis(num).EZ = uValue
        Case "NUXY": Matis(num).NUXY = uValue: Matis(num).Mode2 = 0
        Case "NUYZ": Matis(num).NUYZ = uValue: Matis(num).Mode2 = 0
        Case "NUXZ": Matis(num).NUXZ = uValue: Matis(num).Mode2 = 0
        Case "PRXY": Matis(num).NUXY = uValue: Matis(num).Mode2 = 1
        Case "PRYZ": Matis(num).NUYZ = uValue: Matis(num).Mode2 = 1
        Case "PRXZ": Matis(num).NUXZ = uValue: Matis(num).Mode2 = 1
        Case "GXY": Matis(num).GXY = uValue
        Case "GYZ": Matis(num).GYZ = uValue
        Case "GXZ": Matis(num).GXZ = uValue
        Case Else: nCode = CN_NUL
    End Select
    If Top_Mati < num Then Top_Mati = num
End Function

Public Function Add_Real(num&, R1#)
    Reals(num).R(1) = R1
    If Top_Real < num Then Top_Real = num
End Function

Public Function Mati_Matrix(Mat1 As Ans_Mati)
    Dim H#
    If Mat1.EY = 0 Then '各向同性
        H = Mat1.EX / (1 - Mat1.NUXY ^ 2)
        With Mat1
            .ArrD(1, 1) = H
            .ArrD(1, 2) = H * .NUXY
            .ArrD(1, 3) = 0
            .ArrD(2, 1) = H * .NUXY
            .ArrD(2, 2) = H
            .ArrD(2, 3) = 0
            .ArrD(3, 1) = 0
            .ArrD(3, 2) = 0
            .ArrD(3, 3) = H * (1 - .NUXY) / 2#
        End With
    Else
        
    End If
End Function

Public Function Node_Theta(Point#(), num&)
    Nodes(num).Theta(0) = Point(0)
    Nodes(num).Theta(1) = Point(1)
    Nodes(num).Theta(2) = Point(2)
End Function

Public Function NXYZ(Pr, ByVal I&, ByVal uScl#)
    Pr(CN_X) = NX(I) + uScl * UX(I)
    Pr(CN_Y) = NY(I) + uScl * UY(I)
    Pr(CN_Z) = 0
End Function
Public Function NX#(I&)
    NX = Nodes(I).Point(0)
End Function

Public Function NY#(I&)
    NY = Nodes(I).Point(1)
End Function

Public Function UX#(I&)
    UX = Node2s(I).U(0)
End Function

Public Function UY#(I&)
    UY = Node2s(I).U(1)
End Function

Function MOPER_MULT&(ParR#(), Par1#(), Par2#())
    Rem 矩阵相乘
    Dim M&, N&, P&, I&, J&, K&, Sum#
    M = UBound(Par1, 1): N = UBound(Par1, 2): P = UBound(Par2, 2)
    If M > UBound(ParR, 1) Then M = UBound(ParR, 1)
    If N > UBound(Par2, 1) Then N = UBound(Par2, 1)
    If P > UBound(ParR, 2) Then P = UBound(ParR, 2)
    For I = 1 To M
        For K = 1 To P
            Sum = 0
            For J = 1 To N
                Sum = Sum + Par1(I, J) * Par2(J, K)
            Next
            ParR(I, K) = Sum
        Next
    Next
End Function

Function MOPER_MULT1&(ParR#(), Par1#(), K#)
    Rem 数乘矩阵
    Dim M&, N&, I&, J&
    M = UBound(Par1, 1): N = UBound(Par1, 2)
    For I = 1 To M
        For J = 1 To N
            ParR(I, J) = K * Par1(I, J)
        Next
    Next
End Function

Function MOPER_TURN&(ParR#(), Par1#())
    Rem 矩阵转置
    Dim M&, N&, I&, J&
    M = UBound(Par1, 1): N = UBound(Par1, 2)
    For I = 1 To M
    For J = 1 To N
        ParR(J, I) = Par1(I, J)
    Next: Next
End Function

Function MOPER_SOLV&(ParR#(), Par1#(), Par2#())
    Rem 解线性方程组
    Dim I&, J&, N&, K&
    N = UBound(Par1, 1)
    For I = 1 To N
        ParR(I, 1) = Par2(I, 1)
    Next
    For I = 1 To N
        '第i行每个元素除以A(i,i)
        If Par1(I, I) <> 0 Then
            For K = I + 1 To N
                Par1(I, K) = Par1(I, K) / Par1(I, I)
            Next
            ParR(I, 1) = ParR(I, 1) / Par1(I, I)
            Par1(I, I) = 1
        Else
            '查找以下A(j,i)不为零的行,处理使A(i,i)=1
            For J = I + 1 To N
                If Par1(J, I) <> 0 Then
                    For K = I + 1 To N
                        Par1(I, K) = Par1(I, K) + Par1(J, K) / Par1(J, I)
                    Next
                    ParR(I, 1) = ParR(I, 1) + ParR(J, 1) / Par1(J, I)
                    Par1(I, I) = 1
                    Exit For
                End If
            Next
            If J > N Then
                'XXX MsgArr Par1
                MsgBox "矩阵奇异!", vbExclamation
                MOPER_SOLV = -1
                Exit Function
            End If
        End If
        '第i行以后,每j行:A(j,k)=A(j,k)/A(j,i)-A(i,k)
        For J = I + 1 To N
            If Par1(J, I) <> 0 Then
                For K = I + 1 To N
                    Par1(J, K) = Par1(J, K) / Par1(J, I) - Par1(I, K)
                Next
                ParR(J, 1) = ParR(J, 1) / Par1(J, I) - ParR(I, 1)
                Par1(J, I) = 0
            End If
        Next
    Next
    '最后求解结果
    For I = N To 1 Step -1
        For J = I + 1 To N
            ParR(I, 1) = ParR(I, 1) - Par1(I, J) * ParR(J, 1)
        Next
    Next
    
End Function