上一篇文章中已經講了通過離散求解的流程，目的就是先構造出每個網格的代數方程，代數方程的具體形式如下，碰到的問題就是需要求解三個系數項bc，af，ac

但是這三個系數項就是采用有限體積方法離散獲取的，每個都有自己的特點， 有的需要當前網格的值或者這個值的梯度，有的需要與當前網格相鄰的單元的值。而且前面也列出了這三個公式的具體表達式，如下

前面我們已經知道求解這三項只需要知道gDiff和Tf還有▽?三個就可以了這三個就可以了，而且這三相都是和網格相關的，而且都需要長篇大論的寫，所以我們就一個一個的處理，首先是這個gDiff吧，先給出一個網格的圖片，圖片來自互聯網哈。

一個典型的非結構網格

這里把公式先寫出來

Ef是表面向量的類正交擴散分量，dCF是兩個單元的質心距離，e是兩個網格單元質心向量的單元向量，Sf是面的法向向量,已經知道的是網格中所有的點的坐標都。

首先來完成單元的質心計算吧，當然質心的計算其實也是一個比較漫長的步驟，首先是把每個面的面心計算出來，對于三角形來說每個面的面心直接加權就可以，但是對于多邊形來說還需要轉化為三角形再處理。

1.大致流程是這樣的，通過面上的頂點坐標加和然后除以頂點的數目可以計算到幾何中心

2.將幾何中心與邊上的兩個點相連得到n邊形對應的n個三角形

3.將這n個三角形的面矢量和面積分別計算出來，還有三角形的幾何中心

4.將這n個三角形的幾何中心和對應的面積加權就可以得到面的面心，三角形面積矢量加和除以2就是這個面的面積矢量Sf。

為了省去你的麻煩，我把函數已經定義好了，當然你也可以自己編寫，最好以類的形式寫出來，計算結果也存為類里面的屬性，可以減少計算量。

"""Created on 2022/1/25 10:45@Author: 123@Email: 985040320@qq.com@Describe: **加入文件描述, 要實現的功能, 注意事項等**"""
    def face_center_and_area(vertices):#傳構成這面的點的列表

        # 計算k個點組成的多邊形的幾何中心
        size = len(vertices)
        geo_center = np.sum(vertices, axis=0) / size
        face_total_area = 0.0
        face_normal_vector = np.zeros((3,))
        face_centroid = np.zeros((3,))

        # 形成三角子面.
        # 每一個面都是以邊作為基礎，以面的幾何中心作為頂點然后重建
        for i in range(len(vertices)):
            # 計算子面的點
            p1 = vertices[i]
            p2 = vertices[(i + 1) % size]
            p3 = geo_center

            # 計算子面的幾何中心
            subface_geo_center = np.sum([p1, p2, p3], axis=0) / 3.0

            # 計算子面的面積和面向量，三角形面積等于面矢量的模長除以2
            sf = np.cross((p2 - p1), (p3 - p1))
            area = np.linalg.norm(sf) / 2.0
            # 計算多邊形面的面積和面矢量
            face_normal_vector += sf
            face_total_area += area
            face_centroid = face_centroid + (area * subface_geo_center)

        face_centroid /= face_total_area
        face_normal_vector /= 2.0
        return face_centroid, face_total_area, face_normal_vector

可見通過上面的四個步驟，我們由點的坐標計算出了面的面積，面心，面心矢量。

面解決了，接下來就是計算體心了，類似于多邊形轉化為三角形再進行向量求面積，求面心，求面矢量，多面體需要轉化為像金字塔的多面體（底面的面積面矢量，面心都已知）之后再進行質心的計算

1.首先由上面的面心計算多面體幾何中心，直接加取平均即可

2.通過這個幾何中心和各個面的面心（上一步計算給出的）可以構造一個多面體并計算出他的質心（在基面和幾何中心的0.25處），這里通過面矢量和這兩個點也可以順便把這個有點像金字塔的幾何的體積同時求出。

3.體積加和就是整個網格的體積，第二部求得的質心進行對應的體積加權就可以得到網格的質心。

"""Created on 2022/1/25 10:45@Author: 123@Email: 985040320@qq.com@Describe: **加入文件描述, 要實現的功能, 注意事項等**"""
    def cell_center_and_volume(faces):#傳入面的列表

        # 計算單元的幾何中心，就是簡單的采用所有點的進行平均值計算
        geo_centeroid = np.zeros((3,))
        for face in faces:
            geo_centeroid += face_center_and_area(face)[0]
        geo_centeroid /= len(faces)

        # 根據每個網格面和幾何中心創建類似金字塔的多面體
        element_volume = 0.0
        element_centroid = np.zeros((3,))

        for face in faces:
            face_centroid, face_area, _ = face_center_and_area(face)

            # 四面體的中心落在基面和幾何中心距離的0.25處
            pyramid_centroid = (0.75 * face_centroid) + (0.25 * geo_centeroid)
            pyramid_volume = (1.0 / 3.0) * face_area * np.linalg.norm(face_centroid - geo_centeroid)
            # 借助四面體的中心計算幾何體的體積
            element_volume += pyramid_volume
            # 所有的面計算完成之后用用對應的體積去加一個權
            element_centroid += (pyramid_volume * pyramid_centroid)

        # 最后計算體積加權網格中心
        element_centroid = element_centroid / element_volume
        return element_centroid, element_volume

看現在網格的質心已經求出來了，dCF是兩個單元的質心距離就是這兩個點構成的向量的模嘛，e就是這兩個點構成的向量再除以這個模嘛，sf就是這個面對應的面矢量嘛,所以就可以計算gDiff了。

現在就已經完成了gDiff和Tf還有▽?三個中一個系數的獲取，是不是就是通過網格的拓撲關系獲取的。后面再依次來處理剩余的兩個，同樣還是通過網格的拓撲結構。

喜歡的朋友可以給個關注或者聯系我