0032] (3)從存儲鄰接體單元的集合的哈希表中提取當前結(jié)點包含的體單元，將所有單 元訪問位為true的單元壓入當前單元分組，并設(shè)置這些單元的單元訪問位為true，然后將 所包含單元其它結(jié)點中結(jié)點訪問位為false的結(jié)點依次壓入結(jié)點隊列。
[0033] (4)判斷所有節(jié)點是否完成遍歷，滿足條件轉(zhuǎn)步驟6,否則判斷單元數(shù)是否達到分 組單元數(shù)上限，條件滿足則將當前分組設(shè)置為下一分組。
[0034] (5)從隊列中提取第一個結(jié)點編號，并將該結(jié)點訪問標識設(shè)置為true，轉(zhuǎn)步驟3。
[0035] (6)分組完成后，將各組內(nèi)結(jié)點根據(jù)所屬單元列表是否完全在當前分組內(nèi)進行分 組，滿足條件的為組內(nèi)點，否則為邊界點。對所有結(jié)點依照分區(qū)及優(yōu)先組內(nèi)點，然后編號邊 界點的順序進行重編號。
[0036] (7)根據(jù)新的編號，更新所有的單元結(jié)點編號信息。并將結(jié)點坐標信息和單元結(jié)點 編號信息存儲在一起。
[0037] (8)根據(jù)各單元頂點信息采用等參元法計算得到各單元和標準單元形函數(shù)的變換 關(guān)系矩陣，各結(jié)點質(zhì)量矩陣及單元受力索引。其中等參元法是將各個四面體單元與正四面 體的形函數(shù)之間建立映射關(guān)系，表現(xiàn)形式為一個3*3矩陣，在每一個時間步中計算頂點坐 標發(fā)生變化后，通過該變換關(guān)系獲取新的形函數(shù)矩陣。同時輸入模型參數(shù)楊氏模量、泊松比 及時間步。
[0038] (10)GPU端初始化：將除結(jié)點坐標信息外的其它信息傳送到GPU端。
[0039] (11)根據(jù)手動輸入或是力反饋設(shè)備捕獲的邊界受力變化或是掃描儀或相機捕獲 的可見表面節(jié)點的位移變化，利用Neo-Hookean模型，該模型定義了力與頂點位移之間的關(guān) 系，總體拉格朗日方法，計算每一個時間步的組成四面體單元的各頂點的受力。在多GPU的 并行實現(xiàn)中將四面體單元的頂點坐標信息分成若干段依次傳輸給GPU，同時在每段數(shù)據(jù)傳
輸完成后開啟GPU端該段數(shù)據(jù)的計算函數(shù) 1 '，： 其中S表示第二Piola-Kirchoff應(yīng)力，μ表示剪切模量，J表示變形梯度，I是三階單位 矩陣，L是右柯西格林張量第一不變量，C是右柯西格林張量，δ是材料的體積彈性模量， 通過該公式計算各單元結(jié)點的受力情況。其中分段大小依據(jù)GPU性能動態(tài)設(shè)定。這種分段 處理方式可以減少多個GPU之間的數(shù)據(jù)交互。
[0040] (12)各單元結(jié)點受力計算完成后，通過預(yù)存的單元受力索引計算結(jié)點所受合力。
[0041] (13)將受力信息傳送回CPU端之后，調(diào)用碰撞檢測和響應(yīng)程序進行處理后，采用 中心差分法公式1
各節(jié)點的坐標變化計算出下一 步的結(jié)點位移。其中At表示時間步，Μ表示四面體單元的質(zhì)量，un+1表示下一個時間步的 位移，Rn表示當前時刻的反作用力，？"表示當前時刻所施加的外力，uni表示前一時刻的位 移，un表示當前時刻的位移即步驟二中捕獲的表面節(jié)點的位移變化，更新四面體單元各頂 點的坐標變化。判斷模型是否已經(jīng)收斂穩(wěn)定。條件滿足則中止計算。否則轉(zhuǎn)步驟11繼續(xù) 計算模型變化。
[0042] 如圖3所示，左側(cè)為變形前的模型，右側(cè)為變形后的結(jié)果。
[0043] 本發(fā)明中未詳細闡述的部分屬于本領(lǐng)域技術(shù)人員的公知技術(shù)。
【主權(quán)項】
1. 一種基于多GPU的人腦變形仿真方法，其特征在于包括如下步驟： (1) 初始化，讀入人腦四面體網(wǎng)格模型數(shù)據(jù)，分析四面體網(wǎng)格單元的拓撲結(jié)構(gòu)關(guān)系得到 四面體網(wǎng)格鄰接體單元的集合，依據(jù)所獲得的集合信息和所用GPU核數(shù)目，對讀入的四面 體網(wǎng)格模型原始數(shù)據(jù)進行分組、編號，然后基于等參元法計算各四面體單元的形函數(shù)矩陣， 用于后面在發(fā)生受力或表面位移后內(nèi)部單元的形變計算； (2) 初始化完成后，根據(jù)手動輸入或是力反饋設(shè)備捕獲的邊界受力變化或是掃描儀或 相機捕獲的可見表面節(jié)點的位移變化，利用Neo-Hookean模型，該模型定義了力與頂點位 移之間的關(guān)系，總體拉格朗日方法，計算每一個時間步的組成四面體單元的各頂點的受力， 在實現(xiàn)過程中采用多GPU加速計算過程，進行數(shù)據(jù)分組減少多GPU數(shù)據(jù)交互； (3) 把步驟二中每個GPU核計算出的受力結(jié)果進行合并，并根據(jù)最終計算結(jié)果采用中，其中At表示時間步，M表示四面 體單元的質(zhì)量，un+1表示下一個時間步的位移，Rn表示當前時刻的反作用力，F(xiàn)n表示當前時 刻所施加的外力即步驟二捕獲的表面受力，Unl表示前一時刻的位移，Un表示當前時刻的位 移即步驟二中捕獲的表面節(jié)點的位移變化，更新四面體單元各頂點的坐標變化，并在多CPU 上實現(xiàn)。2. 根據(jù)權(quán)利要求1所述的基于多GPU的人腦變形仿真方法，其特征在于：步驟⑴中 的對讀入的四面體網(wǎng)格模型原始數(shù)據(jù)進行分組、編號，然后基于等參元法計算各四面體單 元的形函數(shù)矩陣，包括步驟如下： (2. 1)遍歷所有四面體網(wǎng)格單元的頂點，為每一個頂點建立所屬單元列表，為每個頂點 建立訪問標志位并設(shè)置為false; (2. 2)輸入分組起始參考頂點，對四面體網(wǎng)格模型單元進行分組，分組單元數(shù)=四面 體總單元數(shù)/所用GPU核數(shù)，首先將起始點的所屬單元加入分組，并將頂點訪問標志設(shè)為 true，判斷當前分組內(nèi)包含四面體是否達到分組單元數(shù)限制，達到則從下一個輸入的四面 體單元頂點開始，將四面體放入下一分組，否則讀取參考點所屬單元列表內(nèi)訪問標志頂點 為false的節(jié)點的所屬單元并加入當前分組中； (2.3) 步驟（2.2)對四面體單元進行分組完成后，將各組內(nèi)頂點根據(jù)所屬單元列表是 否完全在當前分組內(nèi)進行區(qū)分，滿足條件的為組內(nèi)點，否則為邊界點；對所有分組內(nèi)四面體 單元的頂點依照分區(qū)次序從低到高，對這些四面體單元組內(nèi)頂點和邊界點進行重編號； (2.4) 根據(jù)新的編號，更新所有的單元頂點編號信息； (2. 5)根據(jù)各單元頂點信息采用等參元法計算得到各單元和標準單元形函數(shù)的變換關(guān) 系矩陣，等參元法是將各個四面體單元與正四面體的形函數(shù)之間建立映射關(guān)系，表現(xiàn)形式 為一個3*3矩陣，在每一個時間步中計算頂點坐標發(fā)生變化后，通過該變換關(guān)系獲取新的 形函數(shù)矩陣。3. 根據(jù)權(quán)利要求1所述的多GPU的人腦變形仿真方法，其特征在于：步驟⑵中采用 多GPU加速計算過程，進行數(shù)據(jù)分組減少多GPU數(shù)據(jù)交互，包括步驟如下： (3. 1)初始化：初始化完成后，將四面體單元的頂點的坐標信息和單元節(jié)點編號信息 存儲于同一數(shù)據(jù)結(jié)構(gòu)中，保證在單元信息傳送到GPU端后，已傳送四面體單元相關(guān)的頂點 信息也已經(jīng)傳送完畢； (3. 2)流式傳輸：為了充分利用多GPU服務(wù)器上每個GPU核的計算能力及數(shù)據(jù)傳輸帶 寬，把輸入數(shù)據(jù)分成若干組，每次傳輸一組數(shù)據(jù)后，在開啟GPU計算的同時進行下一組數(shù)據(jù) 的傳輸； (3. 3)在GPU端首先使用由Neo-Hookean模型及總體拉格朗日方法推導(dǎo)出的第二 Piola-Kirchoff應(yīng)力其中S表示第二 Piola-KirchofT應(yīng)力，μ表示剪切模量，J表示變形梯度，I是三階單位矩陣，^是右柯西 格林張量第一不變量，C是右柯西格林張量，δ是材料的體積彈性模量，通過該公式計算各 單元的受力情況； (3. 4)通過初始化計算時生成的各四面體單元頂點受力求和求整個四面體網(wǎng)格模型各 頂點的總受力，并將結(jié)果傳送回CPU端。4.根據(jù)權(quán)利要求1所述的多GPU的人腦變形仿真方法，其特征在于：所述步驟⑵中總 體拉格朗日方法實現(xiàn)過程為：方法在預(yù)測模型變形時通過把整個連續(xù)的變形過程離散為若 干個時刻的變形，設(shè)計算開始的時間為t，計算變形的下一個時刻和時間t的時間間隔叫做 時間步，記做At，首先使用頂點的初始坐標、步驟（1)的形函數(shù)矩陣及捕獲的表面位移或 受力計算時刻t+At時每個四面體單元四個頂點的受力情況，然后通過步驟（3)的計算獲 得t+Δt時刻頂點坐標、形函數(shù)矩陣及捕獲的表面位移或受力計算時刻t+2*At時刻t+Δt 時每個四面體單元四個頂點的受力情況以此類推知道兩個時間步之間頂點位移不再發(fā)生 變化。
【專利摘要】本發(fā)明涉及一種基于多GPU的人腦變形仿真方法，使用總體拉格朗日法、中心差分法及非線性有限元Neo-Hookean模型來進行腦部變形的顯式迭代仿真，通過在GPU端的并行化算法實現(xiàn)加速仿真計算的速度，在引入多GPU進一步提高并行能力之后，對原始數(shù)據(jù)采用廣度優(yōu)先搜索方法進行數(shù)據(jù)分組和重新編號最小化各結(jié)點之間的數(shù)據(jù)相關(guān)性，并通過使用額外的數(shù)據(jù)結(jié)構(gòu)使得單個節(jié)點內(nèi)的輸入數(shù)據(jù)及計算過程存在一一對應(yīng)關(guān)系，進而可以使用流式傳輸實現(xiàn)計算和數(shù)據(jù)傳輸?shù)牟⑿谢?，最終利用混合CPU/GPU的多核計算架構(gòu)進行加速，滿足了人腦變形仿真中的精度與速度要求。本發(fā)明充分利用了人腦變形中的細微性以及時空連續(xù)性特點，可以在計算機構(gòu)造的虛擬環(huán)境中逼真地模擬腦部的變形過程。
【IPC分類】G06T17/30, G06T1/20
【公開號】CN105389853
【申請?zhí)枴緾N201510731095
【發(fā)明人】胡勇, 田野, 沈旭昆
【申請人】北京航空航天大學(xué)
【公開日】2016年3月9日
【申請日】2015年11月2日