高手自動略過。

我本人從CAE行業(yè)跳槽到了 EDA行業(yè)，到現(xiàn)在已經(jīng)有兩年了。發(fā)現(xiàn)EDA行業(yè)的跟數(shù)值計算有關(guān)的節(jié)點軟件也采用有限元，有限差分或者邊界元這些個數(shù)值計算方法。于是現(xiàn)在只能再回憶回憶了，有些內(nèi)容都忘記了。

另外特別說明一點，如果不更新的話，我一般不登陸這個賬號。所以留言的朋友我通常沒辦法及時回復(fù)。

還有 python 有限元編程知識你熟悉有限元編程的起點而不是終點。我擱社會上混了得有4年多了，從沒聽說招聘CAE工程師的人只需要會python就行了的。就算你通過python進(jìn)到了公司，大概率后期還得要求你轉(zhuǎn)C/C++ 的。

參考書籍：有限元基礎(chǔ)教程（國際單位制版）第五版

測試案例 （書上例2.2）： 

這些都是比較簡單的 一維單元，所以我省略了很多詳細(xì)有限元編程流程的步驟，僅適用于初學(xué)者。

1. 我這里創(chuàng)建了一個單元類。node_i 是 單元的第一個節(jié)點編號；node_j 是單元的第二個節(jié)點編號；k 是單元的彈簧常數(shù)。

這里的 form_element_k_matrix() 是產(chǎn)生單元系數(shù)矩陣。

2. 根據(jù)題目描述可以看到，這里只有四個單元。

所以我這里：

然后創(chuàng)建了一個單元列表，方便后面管理：

定義邊界條件數(shù)組：

請注意這里的未知量，我全部使用了 -1.0 來代替。

請繼續(xù)往后看，區(qū)分已知節(jié)點和未知節(jié)點的方式將在下面介紹。

根據(jù)我見過的開源軟件內(nèi)部，其區(qū)分是已知節(jié)點還是未知節(jié)點，都是使用的一個數(shù)組進(jìn)行標(biāo)記。

比如我這里采用了 使用 0 標(biāo)記未知節(jié)點，使用 -1 標(biāo)記已知的第一類邊界條件的節(jié)點；使用 -2 標(biāo)記已知的第二類邊界條件，當(dāng)然本案例里頭只有第一類邊界條件，所以這里我不再添加其他標(biāo)記的相應(yīng)處理流程。

這里只有第一個節(jié)點和最后一個節(jié)點是第一類邊界條件，使用 -1 將他們標(biāo)記為已知節(jié)點。

仔細(xì)閱讀過教材的朋友應(yīng)該能知道，整體系數(shù)矩陣中，要去掉已知節(jié)點的行和列上對應(yīng)的元素，將剩下的數(shù)據(jù)組合成新的系數(shù)矩陣。

這里我定義了一個產(chǎn)生整體系數(shù)矩陣的函數(shù)：

3. 如何來將已知節(jié)點對應(yīng)的行和列從系數(shù)矩陣中去掉

當(dāng)前案例的總的節(jié)點個數(shù)是：5； 已知節(jié)點個數(shù)是 2； 那么剩下未知節(jié)點個數(shù)是： 3

我這里采用的方式是這樣：

首先計算出來原系數(shù)矩陣中的待保留數(shù)據(jù)在新系數(shù)矩陣中要移動的行數(shù)。

那么接下來我們需要需要創(chuàng)建方程組的右端項，右端項也就是 力項。待求量是位移項。

那么這里便跟系數(shù)矩陣劃掉的項有關(guān)系了，其實系數(shù)矩陣劃掉項這個流程并不是直接去掉，而是將已知項移動到了方程組的右端。

故這里進(jìn)行移項處理。

移項處理完成之后，抽取出未知節(jié)點上對應(yīng)的數(shù)據(jù)。

那么接下來，關(guān)于 Ax =b 的計算方式，我們產(chǎn)生了 A 和 b。

調(diào)用 scipy 中的模塊即可進(jìn)行求解了。

計算結(jié)果如下：