senne.vb

一个lbm的程序

            Next
        Next
        dFilledMass()
    End Sub

    Private Sub emptyToItf(ByVal a As Int32, ByVal b As Int32)
        Dim k As Byte
        avgrho = 0
        avgu = 0
        avgv = 0
        cnt = 0
        For k = 1 To q
            If flag(a + ei(k, 0), b + ei(k, 1)) >= CONST_INTERFACE Then
                avgrho += rho(a + ei(k, 0), b + ei(k, 1))
                avgu += u(a + ei(k, 0), b + ei(k, 1))
                avgv += v(a + ei(k, 0), b + ei(k, 1))
                cnt += 1
            End If
        Next
        avgrho /= cnt
        avgu /= cnt
        avgv /= cnt
        rho(a, b) = avgrho
        u(a, b) = avgu
        v(a, b) = avgv
        field(a, b) = New InterfaceCell
        field(a, b).f = getfe0(a, b)
        itfArray.Add(New Int32() {a, b})
        tmpArray.Add(New Int32() {a, b})
        mass(a, b) = 0
        flag(a, b) = CONST_TMPS

    End Sub

    Private Sub dFilledMass()
        Dim a, b As Int32
        Dim f(q) As Double
        Dim location(1) As Int32
        Dim itfcell As InterfaceCell
        Dim k As Byte
        For Each location In tmpArray
            flag(location(0), location(1)) = CONST_INTERFACE
        Next
        tmpArray.Clear()
        For Each location In fillArray
            itfcell = field(location(0), location(1))
            vtotal = 0
            For k = 1 To q
                en = itfcell.normalDirection(0) * ei(k, 0) + itfcell.normalDirection(1) * ei(k, 1)
                If en < 0 Then
                    vi(k) = -en '* s(k)
                Else
                    vi(k) = 0
                End If
                vtotal += vi(k)
            Next

            exmass = mass(location(0), location(1)) - rho(location(0), location(1))
            If vtotal <> 0 Then
                For k = 1 To q
                    a = location(0) + ei(k, 0)
                    b = location(1) + ei(k, 1)
                    If flag(a, b) = CONST_INTERFACE Then
                        mass(a, b) = mass(a, b) + exmass * vi(k) / vtotal
                    End If
                Next
            Else
                cnt = getCellNum(location(0), location(1), CONST_INTERFACE)
                For k = 1 To q
                    a = location(0) + ei(k, 0)
                    b = location(1) + ei(k, 1)
                    If flag(a, b) = CONST_INTERFACE Then
                        mass(a, b) += exmass / cnt
                    End If
                Next
            End If
            f = itfcell.f.Clone
            field(location(0), location(1)) = New FluidCell
            field(location(0), location(1)).f = f.Clone
            flag(location(0), location(1)) = CONST_FLUID
            mass(location(0), location(1)) = rho(location(0), location(1))
            eps(location(0), location(1)) = 1
        Next

    End Sub

    Private Sub dempty()
        Dim a, b As Int32
        Dim location(1) As Int32
        Dim icell As InterfaceCell
        Dim f(q) As Double
        Dim k As Byte
        For Each location In emptyArray
            icell = field(location(0), location(1))
            For k = 1 To q
                a = location(0) + ei(k, 0)
                b = location(1) + ei(k, 1)
                If flag(a, b) = CONST_FLUID Then
                    flag(a, b) = CONST_INTERFACE
                    itfArray.Add(New Int32() {a, b})
                    mass(a, b) = rho(a, b)
                    f = field(a, b).f.Clone
                    field(a, b) = New InterfaceCell
                    field(a, b).f = f.Clone
                End If
                en = icell.normalDirection(0) * ei(k, 0) + icell.normalDirection(1) * ei(k, 1)
                If en > 0 Then
                    vi(k) = en '* s(k)
                Else
                    vi(k) = 0
                End If
            Next
            vtotal = 0
            For k = 1 To q
                vtotal += vi(k)
            Next
            If vtotal <> 0 Then
                For k = 1 To q
                    a = location(0) + ei(k, 0)
                    b = location(1) + ei(k, 1)
                    If flag(a, b) = CONST_INTERFACE Then
                        mass(a, b) += mass(location(0), location(1)) * vi(k) / vtotal
                    End If
                Next
            Else
                cnt = getCellNum(location(0), location(1), CONST_INTERFACE)
                For k = 1 To q
                    a = location(0) + ei(k, 0)
                    b = location(1) + ei(k, 1)
                    If flag(a, b) = CONST_INTERFACE Then
                        mass(a, b) += mass(location(0), location(1)) / cnt
                    End If
                Next
            End If
            field(location(0), location(1)) = New BasicCell
            flag(location(0), location(1)) = CONST_EMPTY
            mass(location(0), location(1)) = 0
            rho(location(0), location(1)) = 0
            eps(location(0), location(1)) = 0
            u(location(0), location(1)) = 0
            v(location(0), location(1)) = 0
        Next
    End Sub

    Private Sub deps()
        Dim location(1) As Int32
        For Each location In itfArray
            eps(location(0), location(1)) = mass(location(0), location(1)) / rho(location(0), location(1))
        Next
    End Sub

    Private Sub dcolor()
        Dim location(1) As Int32
        Dim itfcell As InterfaceCell
        For Each location In itfArray
            st = getNormalCellNum(location(0), location(1))
            itfcell = field(location(0), location(1))
            If st(1) = 0 Then
                If mass(location(0), location(1)) < 0.1 * rho(location(0), location(1)) Then
                    itfcell.flag = CONST_EMPTY
                End If
            End If

            If st(0) = 0 Then
                'If mass(x, y) > (0.9 * rho(x, y)) Then
                itfcell.flag = CONST_FLUID
                'End If
            End If
        Next

    End Sub

    Private Sub ddel()
        exmass = 0
        cnt = 0
        Dim location(1) As Int32
        Dim itfcell As InterfaceCell
        Dim f(q) As Double
        Do While cnt < itfArray.Count
            location = itfArray.Item(cnt)
            itfcell = field(location(0), location(1))
            If itfcell.flag <> CONST_INTERFACE Then
                If itfcell.flag = CONST_EMPTY Then
                    exmass += mass(location(0), location(1))
                    mass(location(0), location(1)) = 0
                    eps(location(0), location(1)) = 0
                    u(location(0), location(1)) = 0
                    v(location(0), location(1)) = 0

                    field(location(0), location(1)) = New BasicCell
                Else
                    exmass += (mass(location(0), location(1)) - rho(location(0), location(1)))
                    mass(location(0), location(1)) = rho(location(0), location(1))
                    eps(location(0), location(1)) = 1
                    f = itfcell.f.Clone
                    field(location(0), location(1)) = New FluidCell
                    field(location(0), location(1)).f = f.Clone
                End If
                flag(location(0), location(1)) = itfcell.flag
                itfArray.RemoveAt(cnt)
            Else
                cnt += 1
            End If
        Loop
        vtotal = 0
        If itfArray.Count > 0 Then
            For Each location In itfArray
                vtotal += eps(location(0), location(1))
            Next

            For Each location In itfArray
                mass(location(0), location(1)) += exmass * (eps(location(0), location(1)) / vtotal)
            Next
        End If
    End Sub

End Class