稀疏矩陣處理的案例
基于 MATLAB 的 ANSYS Harwell-Boeing 格式稀疏矩陣提取工具 —— 剛度矩陣與質(zhì)量矩陣 ￥30
在有限元分析中，ANSYS 可以導(dǎo)出大規(guī)模稀疏矩陣（如剛度矩陣、質(zhì)量矩陣），通常使用 Harwell-Boeing (HB) CCS 格式。這些矩陣對(duì)后續(xù)二次開發(fā)、動(dòng)力學(xué)分析或自定義求解器非常重要，但由于其稀疏和壓縮存儲(chǔ)形式，直接在 MATLAB 中讀取和使用并不方便。
本文提供了 兩個(gè) MATLAB 函數(shù)，可直接從 ANSYS 導(dǎo)出的 HB 矩陣文件中讀取并重構(gòu)成 MATLAB 稀疏矩陣：
1.剛度矩陣提取函數(shù)
輸入：ANSYS 導(dǎo)出的剛度矩陣 HB 文件（stiff.txt）
輸出：MATLAB 稀疏矩陣 K，可直接用于動(dòng)力學(xué)計(jì)算或驗(yàn)證
支持自動(dòng)對(duì)稱化，保證數(shù)值正確
2.質(zhì)量矩陣提取函數(shù)
輸入：ANSYS 導(dǎo)出的質(zhì)量矩陣 HB 文件（mass.txt）
輸出：MATLAB 稀疏矩陣 M
使用與剛度矩陣同樣的解析邏輯，無需額外修改
案例說明：
本文以高速鐵路接觸網(wǎng)結(jié)構(gòu)為例，展示了如何將 ANSYS 中導(dǎo)出的稀疏剛度矩陣和質(zhì)量矩陣，在 MATLAB 中完整展開，并進(jìn)行后續(xù)動(dòng)力學(xué)分析準(zhǔn)備。
通過該方法，可以將大規(guī)模有限元矩陣快速轉(zhuǎn)化為 MATLAB 可操作形式，為自定義振動(dòng)分析、模態(tài)分析及其他科研或工程應(yīng)用提供基礎(chǔ)。
優(yōu)勢(shì)與應(yīng)用：
支持大規(guī)模稀疏矩陣解析
自動(dòng)對(duì)稱化，保證數(shù)值精度
適用于剛度矩陣、質(zhì)量矩陣、其他 HB 格式矩陣
可作為動(dòng)力學(xué)求解器或后處理工具的基礎(chǔ)模塊
使用方法：
1.使用以下代碼對(duì)ansys中生成的質(zhì)量及剛度矩陣進(jìn)行提取，file,5,full（5為工作目錄下full文件的文件名，例如：filename.full）。
展開 稀疏矩陣高性能求解器PARDISO的使用
在實(shí)際中，MKL版本的PARDISO性能是非常高的，在某些稀疏矩陣的條件下，求解速度可以遠(yuǎn)超matlab。
以下，是某個(gè)條件數(shù)比較大但是方程組規(guī)模不大（61349個(gè)未知數(shù)）的matlab和PARDISO求解速度對(duì)比：
在該問題中，采用AMD3900X的處理器，MKL_PARDISO求解時(shí)間為43.6s，求解速度是matlab的兩倍多。該稀疏矩陣的下載鏈接如下：
鏈接：https://pan.baidu.com/s/1lw-SgNlnSUAB4o9dwJa14Q
提取碼：xmim
以上，即是高性能稀疏矩陣線性方程組求解器的簡單介紹，感謝您的閱讀，歡迎關(guān)注公眾號(hào) 有限元術(shù)
展開 大規(guī)模稀疏矩陣線性方程組求解可以有多快！
以自由度數(shù)為41w的如下分布的稀疏矩陣為例，
該稀疏矩陣形成的系數(shù)方程組在matlab中采用直接求解為1.5s,采用一定預(yù)處理下的PCG求解時(shí)間為0.55s。
用筆者自行開發(fā)的稀疏矩陣PCG求解器運(yùn)行此算例，方程組求解時(shí)間為0.628s。
此時(shí)，整個(gè)程序的運(yùn)行時(shí)間瓶頸實(shí)際上并不在方程求解，而在于從文件中讀取稀疏矩陣對(duì)應(yīng)的數(shù)據(jù)。一般情況下采用fscanf讀取數(shù)據(jù)會(huì)快于用fstream讀取。改用fscanf讀取數(shù)據(jù)后，程序總運(yùn)行時(shí)間從原來的16s變?yōu)?s。
【完】
歡迎關(guān)注公眾號(hào) 有限元術(shù)
一個(gè)講有限元技術(shù)的公眾號(hào)
展開 從單元連接關(guān)系到節(jié)點(diǎn)鄰接點(diǎn)-有限元采用稀疏矩陣求解的前置工作
在有限元求解中，最終通常要求解的是一個(gè)關(guān)于場(chǎng)變量的線性方程組，在常見的位移場(chǎng)有限元中，要求解的是各個(gè)節(jié)點(diǎn)的位移，該線性方程組的系數(shù)矩陣通常稱為剛度矩陣，方程組右邊通常稱為右端項(xiàng)或者荷載向量。一般情況下，由于網(wǎng)格劃分后并不是所有節(jié)點(diǎn)都兩兩連接，因此實(shí)際上最終形成的整體剛度矩陣中大部分元素為0，這種矩陣稱為稀疏矩陣。在有限元求解中，對(duì)于這種系數(shù)矩陣為稀疏矩陣的方程組，一種常見的方法是僅保存剛度矩陣的非0元素到內(nèi)存中，0元素不保存，這樣就可以以更小的內(nèi)存保存大型結(jié)構(gòu)的剛度矩陣。
那具體矩陣中有多少元素為0，就可以認(rèn)為其是稀疏矩陣呢？這個(gè)界限實(shí)際上比較模糊，有文獻(xiàn)給出如下定義：如果矩陣的A的非0元素?cái)?shù)量為O(n)，其中n是A的階數(shù)，則矩陣為稀疏矩陣。
稀疏矩陣經(jīng)常通過非0元素分布圖表示其稀疏性質(zhì)，以下是兩個(gè)常見的稀疏矩陣的分布圖：
在有限元分析中，非0元素的分布，實(shí)際上主要取決于單元的節(jié)點(diǎn)連接，以下圖中的單元連接為例：
假設(shè)圖中每個(gè)節(jié)點(diǎn)一個(gè)自由度，則整體剛度矩陣為16x16的矩陣，而具體非0元素的分布，可以通過單元連接得到鄰接點(diǎn)得到，所謂鄰接點(diǎn)，指的是相對(duì)于當(dāng)前單元位于同一單元內(nèi)的所有點(diǎn)的集合。以節(jié)點(diǎn)6為例，其鄰接點(diǎn)是1,2,3,5,7,9,10,11。
獲得上述鄰接點(diǎn)后，剛度矩陣中第6行的非0元素的位置實(shí)際上就確定了：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)。
在實(shí)際采用稀疏矩陣求解有限元問題時(shí)，獲得上述非0元素位置后，就可以對(duì)剛度矩陣采用稀疏矩陣存儲(chǔ)，常見的存儲(chǔ)方式有COO，CSR，CSC和DIA等。
展開 
從單元連接關(guān)系到節(jié)點(diǎn)鄰接點(diǎn)-有限元中形成稀疏矩陣求解的前置工作
在有限元求解中，最終通常要求解的是一個(gè)關(guān)于場(chǎng)變量的線性方程組，在常見的位移場(chǎng)有限元中，要求解的是各個(gè)節(jié)點(diǎn)的位移，該線性方程組的系數(shù)矩陣通常稱為剛度矩陣，方程組右邊通常稱為右端項(xiàng)或者荷載向量。一般情況下，由于網(wǎng)格劃分后并不是所有節(jié)點(diǎn)都兩兩連接，因此實(shí)際上最終形成的整體剛度矩陣中大部分元素為0，這種矩陣稱為稀疏矩陣。在有限元求解中，對(duì)于這種系數(shù)矩陣為稀疏矩陣的方程組，一種常見的方法是僅保存剛度矩陣的非0元素到內(nèi)存中，0元素不保存，這樣就可以以更小的內(nèi)存保存大型結(jié)構(gòu)的剛度矩陣。
那具體矩陣中有多少元素為0，就可以認(rèn)為其是稀疏矩陣呢？這個(gè)界限實(shí)際上比較模糊，有文獻(xiàn)給出如下定義：如果矩陣的A的非0元素?cái)?shù)量為O(n)，其中n是A的階數(shù)，則矩陣為稀疏矩陣。
稀疏矩陣經(jīng)常通過非0元素分布圖表示其稀疏性質(zhì)，以下是兩個(gè)常見的稀疏矩陣的分布圖：
在有限元分析中，非0元素的分布，實(shí)際上主要取決于單元的節(jié)點(diǎn)連接，以下圖中的單元連接為例：
假設(shè)圖中每個(gè)節(jié)點(diǎn)一個(gè)自由度，則整體剛度矩陣為16x16的矩陣，而具體非0元素的分布，可以通過單元連接得到鄰接點(diǎn)得到，所謂鄰接點(diǎn)，指的是相對(duì)于當(dāng)前單元位于同一單元內(nèi)的所有點(diǎn)的集合。以節(jié)點(diǎn)6為例，其鄰接點(diǎn)是1,2,3,5,7,9,10,11。
獲得上述鄰接點(diǎn)后，剛度矩陣中第6行的非0元素的位置實(shí)際上就確定了：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)。
在實(shí)際采用稀疏矩陣求解有限元問題時(shí)，獲得上述非0元素位置后，就可以對(duì)剛度矩陣采用稀疏矩陣存儲(chǔ)，常見的存儲(chǔ)方式有COO，CSR，CSC和DIA等。
展開 