基於gpu並行加速的無網格伽遼金法結構拓撲優化方法
2023-05-30 07:00:21 6
基於gpu並行加速的無網格伽遼金法結構拓撲優化方法
【專利摘要】本發明公開了一種基於GPU並行加速的無網格伽遼金法結構拓撲優化方法。本發明主要步驟如下：(1)將數據讀入主機內存，通過CPU布置積分點，建立節點、積分點與局部搜索小格子的關係，計算節點影響域半徑和積分點定義域半徑，確立節點和積分點關係，然後將數據複製到GPU全局存儲器中；(2)根據不同的計算數據設定不同的GPU線程塊與線程數量；(3)通過CPU和GPU異步組裝和求解無網格伽遼金法總體離散系統方程，得到位移近似解；(4)在GPU中進行結構拓撲優化計算，通過設計變量的殘差來判斷迭代是否結束及結果是否輸出。本發明硬體成本低，通用性強，在滿足工程精度要求的前提下可降低大量時間消耗。
【專利說明】基於GPU並行加速的無網格伽遼金法結構拓撲優化方法
【技術領域】
[0001]本發明涉及計算機輔助工程中的結構優化領域,具體涉及一種基於GPU (GraphicsProcessing Unit)並行加速的無網格伽遼金法結構拓撲優化方法。
【背景技術】
[0002]近二十年來，在計算力學領域，無網格法脫穎而出。無網格法不是採用對場變量插值和近似的網格單元，而是採用對場變量進行逼近和近似的節點，同時利用權函數來表徵節點信息。現在已經發展了十餘種無網格方法，如光滑質點流體動力學法、無網格伽遼金法等。無網格伽遼金法是一種非常流行的無網格法，它採用移動最小二乘法構造位移近似形函數，從能量泛函的弱變分形式得到整體控制方程，再利用拉格朗日法、罰函數法等處理本質邊界條件，最後求解總體離散系統方程得到節點位移的數值解。無網格伽遼金法具有較高的協調性和穩定性，計算精度高，收斂快，而且消除了剪切鎖死和體積鎖死等缺點，現已被廣泛應用於衝擊動力學、金屬塑性成型、結構優化等問題，但是其剛度矩陣的組裝和總體離散系統方程求解複雜，且在計算節點影響域半徑、組裝剛度矩陣等時採用全局搜索，導致計算量大，速度慢，耗時長，這些缺點嚴重限制了無網格伽遼金法在結構拓撲優化問題中的應用。
[0003]結構拓撲優化在工程結構設計的初始階段起著重要作用，可為設計者提供一個結構布局最優方案的概念性設計，現已被廣泛應用於建築、機械和航天等領域。結構拓撲優化是使設計域內的材料變為有或者無的過程，即對結構傳遞力作用貢獻小的區域，應刪除材料，而對結構傳遞力作用有益的部位，應保留材料。為求得最佳的材料分配方案，使材料能最大程度的利用，需要構造結構拓撲優化問題的優化模型，確立設計變量，約束條件，以及目標函數，經過計算求得目標函數的最小值。結構拓撲優化可以節約材料，降低成本，但是優化模型求解要經過多次迭代計算，耗時很長。
[0004]隨著計算技術和計算機技術的迅猛發展，並行計算已經在多個領域內取得了成功，成為提高計算效率，節約計算成本的有效方法。並行計算是指在並行計算機上，將一個任務分解成多個子任務，然後分配給不同的處理器，各個處理器之間相互協同，並行地執行子任務，從而達到提高求解速度、增大求解問題規模的目的。現有的並行技術主要有分布式並行計算、並行機、多線程等，主要是以CPU(Central Processing Unit)為計算核心的計算平臺，雖然可以取得較好的加速比，但同時存在以下的不足，首先，以CPU為計算核心的並行機價格昂貴，維護和使用相對複雜，對算法的設計和編程能力要求也較高，較難普及；其次，對於大規模問題，會導致分布式計算機之間或者多線程之間的通訊消耗提升，從而降低了整體加速比。
[0005]2006年，NVIDIA公司公布了統一著色器模型和GeForce8系列GPU，GPU從此進入了通用計算機時代。統一著色器模型整合了頂點著色器和片段著色器，使每個著色器都可以自由承擔原本某種特定的著色器的工作，這種方式使得GPU在圖形處理時空閒的著色器更少，計算效率更高。隨後，NVIDIA公布了統一計算設備架構CUDA(Compute UnifiedDevice Architecture)，並發布了專門針對CUDA架構的程式語言CUDA C。CUDA C支持C語言的大多數指令和語法，並加入了使程序能在GPU上進行多線程計算的擴展語言。NVIDA公司推出的CUDA平臺，作為一種通用的GPU並行計算環境，已成為當今最流行的GPU並行計算平臺。
[0006]隨著GPU的廣泛應用，其強大的浮點運算能力越來越受到人們的重視，基於GPU的通用並行計算技術得到快速發展。在數值代數中，存在較為密集的數據計算，如矢量點積、線性方程組求解等，當處理的規模較大時，計算量是非常可觀的，並且大多數科學與工程計算問題，最後的步驟一般歸結為求解線性方程組，因此研究GPU加速代數計算，尤其是大規模線性方程組的GPU加速求解，其加速比一般可達到15?25倍，但是目前無網格伽遼金法結構拓撲優化仍然不能利用基於CUDA架構的GPU通用計算平臺來實現並行加速。
【發明內容】
[0007]本發明的目的在於針對無網格伽遼金法的結構拓撲優化中存在計算量大、耗時長等問題，提供一種基於GPU並行加速的無網格伽遼金法結構拓撲優化方法。
[0008]本發明基於GPU並行加速的無網格伽遼金法結構拓撲優化方法，包括如下順序的步驟:
[0009](I)將無網格計算模型的節點坐標、有限元背景積分網格數據、局部搜索格子數據、本質邊界數據、懲罰係數、節點影響域係數、材料物理參數、結構拓撲優化參數導入主機內存中；
[0010](2)在CPU中根據主機內存中存儲的模型有限元背景積分網格數據、本質邊界數據，分別為模型內部、力邊界和位移邊界布置積分點，當為二維幾何模型問題時，利用程序自帶的Hammer和Gauss數值積分信息，在模型的有限元三角形背景積分網格和線邊界上布置積分點，當為三維幾何模型問題時，利用程序自帶的Hammer數值積分信息，在模型的有限元四面體背景積分網格和本質邊界上的有限元三角形背景積分網格中布置積分點，求得到積分點的坐標、權係數和雅可比係數；
[0011](3)根據局部搜索格子數據劃分出每個局部搜索小格子的區域，然後分別求得模型節點、積分點和局部搜索小格子的關係數據；
[0012](4)採用局部搜索法計算每個節點的影響域半徑和每個積分點的定義域半徑，然後計算每個節點影響域內的積分點數據，以及每個積分點定義域內的節點數據；
[0013](5)將模型節點坐標、模型內部和位移邊界上的積分點數據、懲罰係數、材料物理參數、節點與局部搜索小格子的關係數據、積分點與局部搜索小格子的關係數據、結構拓撲優化參數、節點影響域內的積分點數據、積分點定義域內的節點數據複製到GPU全局存儲器中；
[0014](6)通過GPU讀取GPU全局存儲器中存儲的模型節點坐標、積分點數據、積分點與局部搜索小格子的關係數據、積分點定義域內的節點數據，然後設置GPU中的線程塊的數量和每一個線程塊所包含的線程數量，並行加速計算每個積分點定義域內各個節點相應的形函數值及形函數導數值；
[0015](7)根據GPU全局存儲器中存儲的積分點數據、積分點定義域內的節點數據、節點影響域內的積分點數據、結構拓撲優化參數、節點形函數值，首先定義模型節點的初始相對密度參數，然後計算結構總體積的初始值及靈敏度；
[0016](8)在GPU執行步驟(6)?(7)時，根據模型節點數據、節點影響域半徑、節點與局部搜索小格子的關係數據、局部搜索格子數據，在CPU中確定影響域有重合區域的兩兩節點，組成交叉節點對，並將交叉節點對信息存儲於數組中；
[0017](9)在CPU中讀取主機內存中存儲的模型節點坐標、本質邊界上的積分點數據、本質邊界數據、懲罰係數，利用罰函數法施加本質邊界條件，在CPU中求得總體外力向量和總體懲罰力向量，然後將交叉節點對數據、總體外力向量、總體懲罰力向量複製到GPU全局存儲器中；
[0018](10)通過GPU進行無網格伽遼金法結構拓撲優化的迭代計算；
[0019](11)判斷是否需要輸出無網格伽遼金法結構拓撲優化的計算結果，如果需要則將GPU全局存儲器中存儲的結構拓撲優化的計算結果複製到主機內存中，再寫入輸出文件。
[0020]具體地說，所述步驟(10)包括如下步驟:
[0021](a)通過GPU讀取GPU全局存儲器中存儲的節點形函數值、模型內部積分點數據、積分點定義域內的節點關係數據、模型內部積分點處的初始相對密度參數，然後建立線程與積分點一一映射模式，並行加速計算模型內部積分點處的相對密度值；
[0022](b)通過GPU讀取GPU全局存儲器中存儲的模型節點坐標、節點的影響域半徑、材料物理參數、積分點數據、節點形函數值及形函數導數值、節點與局部搜索小格子的關係數據、積分點與局部搜索小格子的關係數據、本質邊界數據、懲罰係數，然後建立GPU線程塊與交叉節點對一一映射模式，GPU線程和組成交叉節點對的兩個節點公共影響域內的積分點一一映射模式，並行加速組裝總體剛度矩陣和總體懲罰剛度矩陣，並以按行壓縮存儲CSR格式存儲於GPU全局存儲器中；
[0023](c)通過GPU讀取GPU全局存儲器中存儲的總體剛度矩陣、總體懲罰剛度矩陣、總體外力向量和總體懲罰力向量，調用CUDA架構平臺提供的免費函數庫CUBLAS，採用共軛梯度法求解總體離散系統方程，求得無網格伽遼金法的位移近似解；
[0024](d)通過GPU讀取GPU全局存儲器中存儲的積分點數據、節點影響域內的積分點數據、積分點定義域內的節點數據、節點形函數值及形函數導數值、節點無網格伽遼金法位移近似解、材料物理參數以及結構拓撲優化參數，建立線程塊與節點一一映射模式，線程與節點影響域內的積分點一一映射模式，並行計算目標函數的靈敏度；
[0025](e)利用OC準則，通過GPU讀取GPU全局存儲器中的模型內部積分點處的相對密度值、結構總體積的初始值和靈敏度、目標函數的靈敏度、節點形函數值、積分點數據、積分點定義域內的節點數據，求得相對密度參數和體積比；
[0026](f)計算設計變量前後兩次迭代的殘差，並判斷殘差是否滿足初始殘差限制，如果滿足則迭代結束，如果不滿足，則重新開始步驟(a)。
[0027]本發明與現有技術相比，具有下述優點:
[0028](I)本發明採用三角形積分網格(二維幾何問題)和四面體積分網格(三維幾何問題)作為無網格伽遼金法的背景積分網格，三角形積分網格和四面體積分網格具有強大的適應性，使得本發明方法可廣泛適用於二維、三維任意形狀的幾何模型問題。
[0029](2)本發明利用局部搜索方法，分別建立模型節點、積分點與局部搜索小格子的聯繫，在進行無網格伽遼金法結構拓撲優化計算時，避免了費時的全局搜索，提高了計算效率。
[0030](3)本發明將無網格伽遼金法總體剛度矩陣的組裝和本質邊界條件的施加結合起來，減少了總體懲罰剛度矩陣的組裝時間，提高了計算效率，同時簡化了程序，增加了程序的易讀性和可移植性，另外本發明將總體剛度矩陣和總體懲罰剛度矩陣以按行壓縮存儲CSR格式存儲於GPU全局存儲器中，大大降低了對GPU全局存儲器存儲容量的需求。
[0031](4)本發明將無網格伽遼金法模型節點形函數值及其形函數導數值的計算，總體離散方程的組裝與求解，結構拓撲優化計算等均移植到GPU上進行並行加速計算，而積分點的布置，總體外力向量和總體懲罰力向量計算，數據的讀入和輸出等均在CPU上執行，而GPU和CPU採用異步執行的方式，有效的避免了由於主機內存和GPU全局內存之間數據交換所產生的額外時間。
[0032](5)本發明在採用GPU進行計算時，根據並行計算對像和計算任務的不同，建立不同的線程映射模式，即主要建立三種映射模式:線程與節點 映射，線程與積分點 映射，線程塊與交叉節點對一一映射。
[0033](6)本發明能夠充分利用GPU強大的浮點計算能力，為保證程序的計算效率、簡潔性、可移植性，採用CUDA C編程，同時能夠利用NVIDIA公司提供的免費線性代數數值計算庫CUBLAS，所以本方法計算成本低，通用性強，程序簡潔，具有很好的適應性。
[0034](7)本發明通過建立GPU線程、線程塊的各種映射模式，在保證求解精度的前提下，能夠大幅度節約計算時間，提高計算效率，增大求解規模。利用基於CUDA構架的GPU強大的浮點計算能力，大規模的GPU計算線程，在進行無網格伽遼金法結構拓撲優化計算時一個線程塊或一個線程負責一對交叉節點對、一個節點、一個積分點的相應計算，充分利用了具有上千處理核的GPU並行計算能力。
【專利附圖】
【附圖說明】
[0035]圖1為本發明方法的主流程框圖。
[0036]圖2為本發明方法中計算模型節點與局部搜索小格子的關係數據的流程框圖。
[0037]圖3為本發明方法中計算無網格伽遼金法結構拓撲優化目標函數的靈敏度的流程框圖。
[0038]圖4為本發明實施例的模型示意圖
[0039]圖5為圖4所示實施例的模型節點示意圖。
[0040]圖6為圖4所示實施例的背景積分網格示意圖。
[0041]圖7為圖4所示實施例的結構拓撲優化結果示意圖。
【具體實施方式】
[0042]下面結合附圖和實施例對本發明作進一步詳細的描述。
[0043]參見圖1，本發明基於GPU並行加速的無網格伽遼金法結構拓撲優化方法的實施步驟如下:
[0044](I)將無網格計算模型的節點坐標、有限元背景積分網格數據、局部搜索格子數據、本質邊界數據、懲罰係數、節點影響域係數、材料物理參數、結構拓撲優化參數導入主機內存中；[0045](2)在CPU中根據主機內存中存儲的模型有限元背景積分網格數據、本質邊界數據，分別為模型內部、力邊界和位移邊界布置積分點，當為二維幾何模型問題時，利用程序自帶的Hammer和Gauss數值積分信息，分別在模型的有限元三角形背景積分網格和線邊界上布置積分點，當為三維幾何模型問題時，利用程序自帶的Hammer數值積分信息，在模型的有限元四面體背景積分網格和本質邊界上的有限元三角形背景積分網格中布置積分點，求得到積分點的坐標、權係數及雅可比係數；
[0046](3)參見圖2，根據局部搜索格子數據劃分出每個局部搜索小格子的區域，然後分別求得模型節點、積分點和局部搜索小格子的關係數據；具體的步驟如下:
[0047](a)根據局部搜索格子數據，劃分局部搜索總格子，得到每個局部搜索小格子的範圍，並以數據形式為每個小格子編號；
[0048](b)根據模型節點編號從小到大順序逐個讀取主機內存中存儲的模型節點坐標；
[0049](C)根據模型節點的坐標和局部搜索小格子的範圍，確定該節點所在小格子的編號;
[0050](d)為cellnumnod數組中與該小格子的編號相對應處的元素值加I ;
[0051](e)判斷模型節點是否編歷，如果不是則跳到步驟(b)繼續提取，如果是則繼續下
一步計算；
[0052](f)根據局部搜索小格子的編號和cellnumnod數組中的元素值，計算每個局部搜索小格子內的第一個節點在cellnode數組中的下標值並存儲於cellnode_idx數組中；
[0053](g)判斷是否編曆局部搜索小格子，如果沒有則跳至步驟(f)繼續執行，如果已經編歷，則執行下一步；
[0054](h)根據模型節點編號從小到大順序依次提取節點的坐標和局部搜索小格子的範圍數據，確定該節點所在小格子的編號，並將該節點編號存儲於數據cellnode中的相應位置；
[0055](i)判斷是否編歷模型節點，如果沒有則繼續執行步驟(h)，如果編歷則繼續下一
止/J/ O
[0056](4)採用局部搜索法計算每個節點的影響域半徑和每個積分點的定義域半徑，然後確定每個節點影響域內的積分點數據，以及每個積分點定義域內的節點數據；
[0057](5)將模型節點坐標、模型內部和位移邊界上的積分點數據、懲罰係數、材料物理參數、節點與局部搜索小格子的關係數據、積分點與局部搜索小格子的關係數據、結構拓撲優化參數、節點影響域內的積分點數據、積分點定義域內的節點數據複製到GPU全局存儲器中；
[0058](6)通過GPU讀取GPU全局存儲器中存儲的模型節點坐標、積分點數據、積分點與局部搜索小格子的關係數據、積分點定義域內的節點數據，設置GPU中的線程塊的數量和每一個線程塊所包含的線程數量，並行加速計算每個積分點定義域內各個節點相應的形函數值及形函數導數值；
[0059](7)根據GPU全局存儲器中的積分點數據、積分點定義域內的節點數據、節點影響域內的積分點數據、結構拓撲優化參數、節點形函數值，首先定義模型節點的初始相對密度參數，然後計算結構總體積的初始值和靈敏度；
[0060](8)在GPU執行步驟(6)?(7)時，根據模型節點坐標、節點的影響域半徑、節點與局部搜索小格子的關係數據、局部搜索格子數據，在CPU中確定影響域有重合區域的兩兩結點，組成交叉節點對，並將交叉節點對信息存儲於數組中；
[0061](9)在CPU中讀取主機內存中存儲的模型節點坐標、本質邊界積分點數據、本質邊界數據、懲罰係數，利用罰函數法施加本質邊界條件，在CPU中計算總體外力向量和總體懲罰力向量，然後將交叉節點對數據、總體外力向量、總體懲罰力向量複製到GPU全局存儲器中；
[0062](10)通過GPU進行無網格伽遼金法結構拓撲優化的迭代計算；
[0063](11)判斷是否需要輸出無網格伽遼金法結構拓撲優化的計算結果，如果需要則將GPU全局存儲器上的結構拓撲優化結果複製到主機內存中，然後通過主機內存寫入輸出文件。
[0064]所述步驟(10)的詳細步驟包括:
[0065](a)通過GPU讀取GPU全局存儲器中存儲的節點形函數值、模型內部積分點數據、積分點定義域內的節點數據、模型內部積分點處的初始相對密度參數值，建立線程與積分點--映射模式，並行加速計算模型內部積分點處的相對密度值；
[0066](b)通過GP U讀取GPU全局存儲器中存儲的模型節點坐標、節點的影響域半徑、材料物理參數、積分點數據、節點形函數值及形函數導數值、節點與局部搜索小格子的關係數據、積分點與局部搜索小格子的關係數據、本質邊界數據、懲罰係數，然後建立GPU線程塊與交叉節點對一一映射模式，GPU線程和組成交叉節點對的兩個節點公共影響域內的積分點一一映射模式，加速組裝總體剛度矩陣和總體懲罰剛度矩陣，並以按行壓縮存儲CSR格式存儲於GPU全局存儲器中；
[0067](c)通過GPU讀取GPU全局存儲器中存儲的總體剛度矩陣、總體懲罰剛度矩陣、總體外力向量、總體懲罰力向量，調用CUDA平臺提供的免費函數庫GUBLAS，採用共軛梯度法求解總體離散系統方程，得到模型各個節點的無網格伽遼金法位移近似解；
[0068](d)參見圖3，通過GPU讀取GPU全局存儲器中存儲的積分點數據、節點影響域內的積分點數據、積分點定義域內的節點數據、節點形函數值及形函數導數值、節點無網格伽遼金法的位移近似解、材料物理參數以及結構拓撲優化參數，建立線程塊與節點一一映射模式，線程與節點影響域內的積分點一一映射模式，並行加速計算目標函數的靈敏度；具體步驟如下:
[0069]第一步:通過GPU從按小到大逐個提取GPU全局存儲器中存儲的模型節點編號，根據節點編號和節點影響域內積分點數據得到該節點影響域內的積分點數和積分點編號；
[0070]第二步:依次提取該節點影響域內的積分點，根據積分點編號和積分點定義域內的節點數據，計算得到該節點影響域內的積分的定義域內的節點編號；
[0071]第三步:通過GPU讀取GPU全局存儲器中存儲的節點形函數值、懲罰因子、內部積分點處的相對密度值、材料物理參數、積分點的權係數和雅可比係數、積分點的編號，計算該積分點處相對密度的(懲罰因子-1)次冪；
[0072]第四步:編歷該節點影響域內的積分點定義域內的節點，根據該節點的編號和GPU全局存儲器中存儲的節點形函數導數值、節點位移近似值、該積分點處相對密度的(懲罰因子-1)次冪值，計算靈敏度分量值；
[0073]第五步:判斷是否編歷該節點影響域內的積分點定義域內的節點，如果沒有，則跳至第四步繼續執行，如果編歷，則執行下一步；
[0074]第六步:累加靈敏度分量，得到靈敏度值；
[0075]第七步:判斷是否編歷該模型節點內的積分點，如果沒有，則跳至第二步繼續執行，如果編歷則執行下一步；
[0076]第八步:判斷是否編歷模型節點，如果沒有則跳至第一步繼續執行，如果編歷則執行下一步。
[0077](e)利用OC (Optimization Criteria)準則,通過GPU讀取GPU全局存儲器中存儲的模型內部積分點處的相對密度值、結構總體積的初始值和靈敏度、目標函數的靈敏度、節點形函數值、積分點數據、積分點定義域內的節點數據，求得相對密度參數和體積比；
[0078](f)計算設計變量的前後兩次迭代的殘差，並判斷殘差是否滿足初始殘差限制，如果滿足則迭代結束，如果不滿足，則重新開始步驟(a)。
[0079]下面是將本發明方法應用於工程上的一個實例:
[0080]參見圖4，本實施例為長15m，寬10m，厚度為lm，受集中力P = IN作用的懸臂梁，材料的彈性E = L 0，泊松比μ = 0.3，無網格伽遼金法結構拓撲優化模型包括1033個節點，其分布如圖5所示，有限元背景積分網格包括1033個節點和1944個三角形網格，如圖6所示。本發明方法針對本實施例的具體實施步驟如下:
[0081]第一步:將已經準備好的模型數據(包括模型的節點坐標、有限元背景積分網格頂點的編號和坐標、本質邊界節點對中兩節點的編號、本質邊界信息、懲罰係數)、材料物理參數(包括材料彈性模量、泊松比)、節點影響域係數、局部搜索格子數據、拓撲優化參數(懲罰因子、最大體積保留率)輸入到input_para.txt文件中，在程序運行時，系統自動將以上數據讀入到主機內存中。一個標準的input_para.txt文件格式如下表一:
[0082]表一
[0083]產模型節點坐標*/
0.0-5.0
產有限元背景積分網格頂點的編號*/
9651029M(
嚴有限元背景積分網格頂點的坐標*/
0.0-5.0
產位移邊界節點對中兩節點的編號*/
6298
產集中力受力節點的編號iV2
產集中力大小*/
I
產材料彈性模量*/
1.0
產材料泊松比*/
0.3
產節點的影響域係數*/
3.0
/*懲罰系數*/
1.0E4
/*局部搜索總格子的邊長1510
，局部搜索總格子X h N的等分數*/』
31
產局部搜索總格子的屮心點坐標*/
7.5O
/*懲罰因子*/.10
產最大體積保留率*/
0.3
[0084]第二步:依據主機內存中存儲的局部搜索總格子邊長和X軸方向的等分數，計算每個局部搜索小格子的邊長；本實施例中局部搜索總格子X軸方向的邊長為15m，y軸方向的邊長為10m，X軸方向的等分數為30，所以本實施例中每個小格子X軸方向的邊長為
0.5m, y軸方向的邊長可根據局部搜索總格子邊長的等分數與邊長的長度成比例求得約為
0.33333m。
[0085]第三步:在CPU中根據本實施例三角形有限元背景積分網格頂點的坐標與編號、程序自帶的Hammer數值積分數據，在模型內部布置積分點，求得積分點的坐標、權係數及雅可比係數，用於存儲模型內部積分點信息的數組需要的存儲空間大小為:1944(本實施例中有限元背景積分網格個數)*3(每個三角形有限元背景積分網格內布置的積分點數)*4*sizeof (double)字節。
[0086]第四步:通過CPU讀取主機內存中存儲的有限元背景積分網格的頂點坐標與編號，本質邊界節點對中兩節點的編號，根據程序自帶的Gauss積分數據，在本質邊界上布置積分點，求得位移和力邊界上積分點的坐標、權係數、雅克比係數，現在分配給位移邊界積分點數組的存儲空間為:24 (本實施例中位移邊界節點對個數)*4 (每個單元節點對上布置的積分點數)*4*sizeof (double)字節，因為本實施例中施加的是集中力，只在受力點布置一個積分點，所以為力邊界積分點數組分配的存儲空間大小為:4*siZe0f (double)字節。
[0087]第五步:如圖2所示，在CPU中，根據模型節點坐標、局部搜索格子數據，建立節點和局部搜索小格子的對應關係，並將其關係數據存儲於cellnode, cellnumnod, cellnode_idx三個數組中，同理建立模型積分點與局部搜索小格子的對應關係。
[0088]第六步:在CPU中，對模型節點進行循環求得節點的動態影響域半徑。根據節點坐標，確定該節點鄰近的幾個局部搜索小格子在全部格子中的編號，並根據局部搜索小格子與節點的對應關係數據確定該節點幾個鄰近局部搜索小格子中所有節點的編號，再在幾個鄰近小格子中進行全局搜索，求得鄰近局部搜索小格子中其它節點和該節點的最小距離rmin，然後利用公式:節點影響域半徑=3.0(節點影響域係數)*rmin可求得該節點的動態影響域半徑，在主機內存中為其分配的存儲空間為:1033 (模型節點的個數)*2 (x、y軸兩個方向)*sizeof (double)字節，同理求得積分點的定義域半徑。
[0089]第七步:在CPU中，讀取模型的節點坐標、積分點坐標、節點與局部搜索小格子的關係數據、積分點與局部搜索小格子的關係數據、節點的影響域半徑，採用局部搜索方法求得每個節點影響域內的積分點數據，每個積分點定義域內的節點數據。 [0090]第八步:將模型節點坐標、模型內部和位移邊界的積分點數據、位移邊界數據、節點與局部搜索小格子的關係數據、積分點與局部搜索小格子的關係數據、結構拓撲優化參數、節點影響域內的積分點數據、積分點定義域內的節點數據複製到GPU全局存儲器中，並為它們分配和主機內存中同樣大小的存儲空間。
[0091]第九步:通過GPU讀取GPU全局存儲器中存儲的模型節點坐標、積分點數據、積分點與局部搜索小格子的關係數據、積分點定義域內的節點數據，建立線程與積分點一一映射模式，通過最小二乘逼近，並行求得每個積分點定義域內節點的形函數值，以及形函數對應於各個坐標軸的導數值，在本實施例中為形函數對x、y軸的導數值。
[0092]第十步:根據GPU全局存儲器中存儲的積分點數據、積分點定義域內的節點數據、節點影響域內的積分點數據、結構拓撲優化參數、節點形函數值，首先建立線程和節點一一映射模式，定義節點的初始相對密度參數，然後建立線程和積分點一一映射模式，算得結構總體積的初始值和靈敏度。
[0093]第H^一步:GPU執行步驟九、十時，在CPU中根據節點坐標、節點的影響域半徑、節點與局部搜索小格子的關係數據、局部搜索格子數據，採用局部搜索法，在CPU中確定影響域有重合區域的兩兩節點，組成交叉節點對，然後將交叉節點對信息以數據形式存儲於主機內存中。
[0094]第十二步:讀取主機內存中存儲的模型節點坐標、本質邊界數據、有限元背景積分網格數據、力邊界積分點信息、節點與局部搜索小格子的關係數據、積分點與局部搜索小格子的關係數據，採用局部搜索法，求得模型的總體外力向量F。
[0095]第十三步:在CPU中讀取主機內存中存儲的模型節點坐標、位移邊界積分點數據、本質邊界數據、懲罰係數、節點與局部搜索小格子的關係數據、積分點與局部搜索小格子的關係數據，利用罰函數法施加本質邊界條件，求得總體懲罰力向量Fa，然後將交叉節點對信息、總體外力向量和總體懲罰力向量複製到GPU全局存儲器中。
[0096]第十四步:讀取GPU全局存儲器中模型節點的形函數值及形函數導數值、積分點數據、材料相對密度值，利用最小二乘逼近，在GPU中建立線程與積分點一一映射模式，計算模型內部積分點處的相對密度值。
[0097]第十五步:在GPU中建立線程塊和交叉節點對一一映射模式，線程和組成交叉節點對的兩個節點公共影響域內的積分點一一映射模式，並行加速組裝總體剛度矩陣K和總體懲罰剛度矩陣κ%並以CSR格式存儲於GPU全局存儲器中，然後讀取GPU全局存儲器中存儲的總體外力向量F和總體懲罰力向量F%在GPU中調用CUDA平臺提供的免費函數，採用共軛梯度法求解控制方程(K+Ka)U = F+Fa，求得無網格伽遼金法的位移近似解U。
[0098]第十六步:如圖3所示，在GPU中計算目標函數的靈敏度值。首先建立線程塊與節點一一映射模式，求得每個節點影響域內的積分點數以及積分點編號，然後建立線程塊內的線程與該節點影響域內的積分點一一映射模式，求得目標函數的靈敏度值。
[0099]第十七步:讀取GPU全局存儲器中存儲的模型節點初始相對密度參數、結構總體積的初始值和靈敏度值、目標函數的靈敏度值、節點的形函數值、積分點數據、積分點定義域內的節點數據，建立積分點與GPU線程一一映射模式，在GPU上利用OC準則計算節點的相對密度值。
[0100]第十八步:在GPU上求得模型節點相對密度的殘差，並判斷殘差是否小於初始設定值，如果大於則跳到第十四步繼續迭代，如果小於則跳出迭代計算，並將GPU全局存儲器中存儲的無網格伽遼金法結構拓撲優化的計算結果複製到主機內存中，再寫入輸出文件，本實施例無網格伽遼金法結構拓撲優化結果如圖7所示。
[0101]雖然參考優先實施例對本發明進行描述，但以上所述並不限定本發明的保護範圍，凡屬於本發明思路下的技術方案均屬於本發明的保護範圍，任何在本發明的精神及原則內的修改、改進等，也應視為本發明的保護範圍。
【權利要求】
1.一種基於GPU並行加速的無網格伽遼金法結構拓撲優化方法，其特徵在於包括如下順序的步驟: (1)將無網格計算模型的節點坐標、有限元背景積分網格數據、局部搜索格子數據、本質邊界數據、懲罰係數、節點影響域係數、材料物理參數、結構拓撲優化參數導入主機內存中； (2)在CPU中根據主機內存中存儲的模型有限元背景積分網格數據、本質邊界數據，分別為模型內部、力邊界和位移邊界布置積分點，當為二維幾何模型問題時，利用程序自帶的Hammer和Gauss數值積分信息，在模型的有限元三角形背景積分網格和線邊界上布置積分點，當為三維幾何模型問題時，利用程序自帶的Hammer數值積分信息，在模型的有限元四面體背景積分網格和本質邊界上的有限元三角形背景積分網格中布置積分點，求得到積分點的坐標、權係數和雅可比係數； (3)根據局部搜索格子數據劃分出每個局部搜索小格子的區域，然後分別求得模型節點、積分點和局部搜索小格子的關係數據； (4)採用局部搜索法計算每個節點的影響域半徑和每個積分點的定義域半徑，然後計算每個節點影響域內的積分點數據，以及每個積分點定義域內的節點數據； (5)將模型節點坐標、模型內部和位移邊界上的積分點數據、懲罰係數、材料物理參數、節點與局部搜索小格子的關係數據、積分點與局部搜索小格子的關係數據、結構拓撲優化參數、節點影響域內的積分點數據、積分點定義域內的節點數據複製到GPU全局存儲器中； (6)通過GPU讀取GPU全局存儲器中存儲的模型節點坐標、積分點數據、積分點與局部搜索小格子的關係數據、積分點定義域內的節點數據，然後設置GPU中的線程塊的數量和每一個線程塊所包含的線程數量，並行加速計算每個積分點定義域內各個節點相應的形函數值及形函數導數值； (7)根據GPU全局存儲器中存儲的積分點數據、積分點定義域內的節點數據、節點影響域內的積分點數據、結構拓撲優化參數、節點形函數值，首先定義模型節點的初始相對密度參數，然後計算結構總體積的初始值及靈敏度； (8)在GPU執行步驟(6)~(7)時，根據主機內存中存儲的模型節點坐標、節點影響域半徑、節點與局部搜索小格子的關係數據、局部搜索格子數據，在CPU中確定影響域有重合區域的兩兩節點，組成交叉節點對，並將交叉節點對信息存儲於數組中； (9)在CPU中讀取主機內存中存儲的模型節點坐標、本質邊界上的積分點數據、本質邊界數據、懲罰係數，利用罰函數法施加本質邊界條件，在CPU中求得總體外力向量和總體懲罰力向量，然後將交叉節點對數據、總體外力向量、總體懲罰力向量複製到GPU全局存儲器中； (10)通過GPU進行無網格伽遼金法結構拓撲優化的迭代計算； (11)判斷是否需要輸出無網格伽遼金法結構拓撲優化的計算結果，如果需要則將GPU全局存儲器中存儲的結構拓撲優化的計算結果複製到主機內存中，再寫入輸出文件。
2.根據權利要求1所述基於GPU並行加速的無網格伽遼金法結構拓撲優化方法，其特徵在於所述步驟(10)包括如下步驟: (a)通過GPU讀取GPU全局存儲器中存儲的節點形函數值、模型內部積分點數據、積分點定義域內的節點關係數據、模型內部積分點處的初始相對密度參數，然後建立線程與積分點—映射模式，並行加速計算模型內部積分點處的相對密度值； (b)通過GPU讀取GPU全局存儲器中存儲的模型節點坐標、節點的影響域半徑、材料物理參數、積分點數據、節點形函數值及形函數導數值、節點與局部搜索小格子的關係數據、積分點與局部搜索小格子的關係數據、本質邊界數據、懲罰係數，然後建立GPU線程塊與交叉節點對一一映射模式，GPU線程和組成交叉節點對的兩個節點公共影響域內的積分點一一映射模式，並行加速組裝總體剛度矩陣和總體懲罰剛度矩陣，並以按行壓縮存儲CSR格式存儲於GPU全局存儲器中； (c)通過GPU讀取GPU全局存儲器中存儲的總體剛度矩陣、總體懲罰剛度矩陣、總體外力向量和總體懲罰力向量，調用CUDA架構平臺提供的免費函數庫CUBLAS，採用共軛梯度法求解總體離散系統方程，求得無網格伽遼金法的位移近似解； (d)通過GPU讀取GPU全局存儲器中存儲的積分點數據、節點影響域內的積分點數據、積分點定義域內的節點數據、節點形函數值及形函數導數值、節點無網格伽遼金法位移近似解、材料物理參數以及結構拓撲優化參數，建立線程塊與節點一一映射模式，線程與節點影響域內的積分點一一映射模式，並行計算目標函數的靈敏度； (e)利用OC準則，通過GPU讀取GPU全局存儲器中的模型內部積分點處的相對密度值、結構總體積的初始值和靈敏度、目標函數的靈敏度、節點形函數值、積分點數據、積分點定義域內的節點數據，求得相對密度參數和體積比； (f)計算設計變量前後兩次迭代的殘差，並判斷殘差是否滿足初始殘差限制，如果滿足則迭代結束，如果不 滿足，則重新開始步驟(a)。
3.根據權利要求2所述基於GPU並行加速的無網格伽遼金法結構拓撲優化方法，其特徵在於所述步驟(d)包括如下步驟: 第一步:通過GPU從按小到大逐個提取GPU全局存儲器中存儲的模型節點編號，根據節點編號和節點影響域內積分點數據得到該節點影響域內的積分點數和積分點編號； 第二步:依次提取該節點影響域內的積分點，根據積分點編號和積分點定義域內的節點數據，計算得到該節點影響域內的積分的定義域內的節點編號； 第三步:通過GPU讀取GPU全局存儲器中存儲的節點形函數值、懲罰因子、內部積分點處的相對密度值、材料物理參數、積分點的權係數和雅可比係數、積分點的編號，計算該積分點處相對密度的(懲罰因子-1)次冪； 第四步:編歷該節點影響域內的積分點定義域內的節點，根據該節點的編號和GPU全局存儲器中存儲的節點形函數導數值、節點位移近似值、該積分點處相對密度的(懲罰因子-1)次冪值，計算靈敏度分量值； 第五步:判斷是否編歷該節點影響域內的積分點定義域內的節點，如果沒有，則跳至第四步繼續執行，如果編歷，則執行下一步； 第六步:累加靈敏度分量，得到靈敏度值； 第七步:判斷是否編歷該模型節點內的積分點，如果沒有，則跳至第二步繼續執行，如果編歷則執行下一步； 第八步:判斷是否編歷模型節點，如果沒有則跳至第一步繼續執行，如果編歷則執行下一步。
4.根據權利要求1所述基於GPU並行加速的無網格伽遼金法結構拓撲優化方法，其特徵在於所述步驟(3)包括如下步驟: (a)根據局部搜索格子數據，劃分局部搜索總格子，得到每個局部搜索小格子的範圍，並以數據形式為每個小格子編號； (b)根據模型節點編號從小到大順序逐個讀取主機內存中存儲的模型節點坐標； (C)根據模型節點的坐標和局部搜索小格子的範圍，確定該節點所在小格子的編號； (d)為cellnumnod數組中與該小格子的編號相對應處的元素值加I; (e)判斷模型節點是否編歷，如果不是則跳到步驟(b)繼續提取，如果是則繼續下一步計算; (f)根據局部搜索小格子的編號和cellnumnod數組中的元素值，計算每個局部搜索小格子內的第一個節點在cellnode數組中的下標值並存儲於cellnode_idx數組中； (g)判斷是 否編曆局部搜索小格子，如果沒有則跳至步驟(f)繼續執行，如果已經編歷，則執行下一步； (h)根據模型節點編號從小到大順序依次提取節點的坐標和局部搜索小格子的範圍數據，確定該節點所在小格子的編號，並將該節點編號存儲於數據cellnode中的相應位置； (i)判斷是否編歷模型節點，如果沒有則繼續執行步驟(h)，如果編歷則繼續下一步。
【文檔編號】G06F17/50GK103970960SQ201410221811
【公開日】2014年8月6日 申請日期:2014年5月23日 優先權日:2014年5月23日
【發明者】龔曙光, 劉奇良, 盧海山, 唐芳, 王恆宇, 王龍 申請人:湘潭大學