在有限元求解中，最終通常要求解的是一個關(guān)于場變量的線性方程組，在常見的位移場有限元中，要求解的是各個節(jié)點的位移，該線性方程組的系數(shù)矩陣通常稱為剛度矩陣，方程組右邊通常稱為右端項或者荷載向量。一般情況下，由于網(wǎng)格劃分后并不是所有節(jié)點都兩兩連接，因此實際上最終形成的整體剛度矩陣中大部分元素為0，這種矩陣稱為稀疏矩陣。在有限元求解中，對于這種系數(shù)矩陣為稀疏矩陣的方程組，一種常見的方法是僅保存剛度矩陣的非0元素到內(nèi)存中，0元素不保存，這樣就可以以更小的內(nèi)存保存大型結(jié)構(gòu)的剛度矩陣。

那具體矩陣中有多少元素為0，就可以認為其是稀疏矩陣呢？這個界限實際上比較模糊，有文獻給出如下定義：如果矩陣的A的非0元素數(shù)量為O(n)，其中n是A的階數(shù)，則矩陣為稀疏矩陣。

稀疏矩陣經(jīng)常通過非0元素分布圖表示其稀疏性質(zhì)，以下是兩個常見的稀疏矩陣的分布圖：

在有限元分析中，非0元素的分布，實際上主要取決于單元的節(jié)點連接，以下圖中的單元連接為例：

假設(shè)圖中每個節(jié)點一個自由度，則整體剛度矩陣為16x16的矩陣，而具體非0元素的分布，可以通過單元連接得到鄰接點得到，所謂鄰接點，指的是相對于當前單元位于同一單元內(nèi)的所有點的集合。以節(jié)點6為例，其鄰接點是1,2,3,5,7,9,10,11。

獲得上述鄰接點后，剛度矩陣中第6行的非0元素的位置實際上就確定了：k(6,1),k(6,2),k(6,3),k(6,5),k(6,6),k(6,7),k(6,9),k(6,10),k(6,11)。

在實際采用稀疏矩陣求解有限元問題時，獲得上述非0元素位置后，就可以對剛度矩陣采用稀疏矩陣存儲，常見的存儲方式有COO，CSR，CSC和DIA等。上面的網(wǎng)格數(shù)量較少，因此可以通過觀察獲得節(jié)點的鄰接點，在實際有限元求解中，網(wǎng)格劃分以后通常得到的是單元的節(jié)點連接信息，即各個單元分別有哪幾個節(jié)點構(gòu)成，并不能直接獲得各個節(jié)點的鄰接點，在上圖中，網(wǎng)格劃分后得到的單元的節(jié)點連接信息如下：

1    1    2    6    5

2    2    3    7    6

3    3    4    8    7

4    5    6   10    9

5    6    7   11   10

6    7    8   12   11

7    9   10   14   13

8   10   11   15   14

9   11   12   16   15

通過上述連接信息，就可以得到任一節(jié)點的鄰接點，在這里，提供一個從上述節(jié)點連接獲得節(jié)點鄰接點的代碼，節(jié)點連接信息按照上面的格式保存在element.txt中。

#include<iostream>
#include<vector>
#include <algorithm>
#include <fstream>
#include<map>
int main(int argc, char** argv)
{
     int nNode ;
     std::cout << "節(jié)點總數(shù)：" << std::endl;
     std::cin >> nNode;
    std::multimap<int,int> *nodeNear=new std::multimap<int,int>[nNode];
    int nElem;
    std::cout << "單元總數(shù)：" << std::endl;
    std::cin >> nElem;
    int elemNumber;
    std::vector<int> *elemConnect=new std::vector<int>[nElem];
    std::ifstream eFile("element.txt");
    if (!eFile) std::cerr << "open file failure" << std::endl;
    
    const int nNodePe=4;   
    int nodeN[nNodePe+1];
    for (int i = 0; i < nElem; i++)
    {
        for (int j = 0; j < nNodePe+1; j++)
        {
            eFile >> nodeN[j];
            elemConnect[i].push_back(nodeN[j]);
        }
        
    }
    for (int i = 0; i <nNode; ++i)
    {
        for (int j = 0; j<nElem; ++j)
        {
            auto a1 = std::find(elemConnect[j].begin()+1, elemConnect[j].end(), i);
            if (a1 != elemConnect[j].end())
            {
                for (auto it = elemConnect[j].begin()+1; it!= elemConnect[j].end(); it++)
                {
                    nodeNear[i].insert(std::pair<int,int>(*it,elemConnect[j][0]));
                }    
            }
        }    
    }
    std::cout << "輸入節(jié)點號：" << std::endl;
    int nodeNumber;
    std::cin >> nodeNumber;
    std::cout << "鄰近節(jié)點：" << std::endl;
    for (auto it:nodeNear[nodeNumber])
    {
        std::cout << "node  " << it.first <<"  in element  "<< it.second << std::endl;
    }
    delete[] nodeNear; delete[] elemConnect;
    return 0;
}

運行結(jié)果如下：

以上，就是從單元連接關(guān)系到節(jié)點鄰接點的具體過程和代碼，感謝您的閱讀！

【完】

歡迎關(guān)注公眾號  有限元術(shù)