[發明專利]基于完全矩陣計算的離散元鄰居搜索及求解方法和系統

在審

• 申請號: 201911066584.1
• 申請日: 2019-11-04
• 公開（公告）號: CN110781448A
• 公開（公告）日: 2020-02-11
• 發明（設計）人: 劉春;劉輝;施斌;王寶軍;樂天呈;張丹;唐朝生;朱鴻鵠;顧凱
• 申請（專利權）人: 南京大學
• 主分類號: G06F17/16
• 分類號: G06F17/16;G06F16/901
• 代理公司: 32249 南京瑞弘專利商標事務所(普通合伙)
• 代理人: 彭雄

1.一種基于完全矩陣計算的離散元鄰居搜索方法，其特征在于，包括以下步驟：

步驟1，單元編號及信息存儲：將單元順序編號，記從1到m，同時添加虛單元，其編號m+1，通過編號唯一標識單元；通過虛單元填充得到接觸搜索的鄰居矩陣表示，鄰居矩陣為m行，每一行記錄單元可能接觸的單元編號，長短不一部分用虛單元編號填充，單元信息存儲為單元屬性數組，記為aX、aY、aZ和aR，其中，編號m+1單元為虛單元，aX、aY、aZ為三維坐標，aR為半徑；

步驟2，網格劃分及編號：

步驟21，輸入不包括虛單元的單元三維坐標及半徑數組，記三維坐標數組為X,Y,Z，半徑數組為R，其均為m行的一維數組；

步驟22，根據三維坐標數組確定網格邊界及邊長：

p方向網格兩側邊界：

Pmin＝min(P)-minR/2

Pmax＝max(P)+minR/2

其中，[Pmin，Pmax]為p方向兩側邊界，[min(P),max(P)]為p方向坐標數組最小最大值，minR為最小單元半徑，P表示x,y,z方向坐標數組；

記[Xmin,Xmax],[Ymin,Ymax],[Zmin,Zmax]為x,y,z方向網格兩側邊界；

網格元胞邊長：

gSide＝2*maxR+dSide

其中，gSide為網格元胞邊長，maxR為單元最大半徑，dSide為可調值；

以最小值邊界為起點，gSide為網格元胞邊長，劃分空間網格；則p方向含虛元胞的元胞數PNum：

其中，ceil表示向上取整運算，元胞沿p方向順序編號，記從1到PNum；

記x,y,z方向元胞數分別為xNum,yNum,zNum,則總元胞數cellNum：

cellNum＝xNum*yNum*zNum

步驟23，計算單元-元胞三維下標索引：

單元所在元胞在p方向的下標索引：

其中，floor表示向下取整運算；

記單元所在元胞沿x,y,z方向所在元胞索引為xIndex,yIndex,zIndex；

步驟3，坐標映射：將元胞沿x->y->z的順序編號，記為1到cellNum，則由步驟23單元-元胞三維下標索引得單元-元胞列表：

aIndex＝xIndex+xNum*(yIndex-1+yNum*(zIndex-1))

其中，aIndex表示單元-元胞列表；

步驟4，索引轉換：根據步驟23得到的單元-元胞三維下標索引，得到鄰域元胞三維下標索引：

[dcxIndex,dcyIndex,dczIndex]＝ndgrid([-1,0,1],[-1,0,1],[-1,0,1])

cxIndex＝xIndex’+dcxIndex(:)

其中，dcxIndex,dcyIndex,dczIndex為元胞三維鄰域算子，大小為3×3×3的數組，其元素取值為-1，0或1，ndgrid為生成格點運算，xIndex’表示xIndex的轉置，dcxIndex(:)表示將其中元素排成一列；記單元鄰域元胞三維下標索引為cxIndex,cyIndex,czIndex；

則由索引轉換，單元-鄰域元胞矩陣：

cIndex＝cxIndex+xNum*(cyIndex-1+yNum*(czIndex-1))

其中，cIndex表示單元-鄰域元胞矩陣；

步驟5，建立元胞-包含單元矩陣：

步驟51，根據步驟3得到的單元-元胞列表，按元胞排序得到元胞-單元列表：

[sortaIndex,sortI]＝sort(aIndex)

其中，sort為排序運算，sortaIndex為排序后的元胞列表，sortI為排序后的對應單元索引；

步驟52，根據排序后的元胞-單元列表生成行索引，記為rowIndex；

步驟53，索引轉換及索引賦值：目標矩陣元胞-包含單元矩陣為一m列的矩陣，其每一列存儲該元胞所包含的單元編號，對于虛元胞及未填滿的列，可用編號0填充；元胞-包含單元矩陣初始化：

cellBallList＝zeros(maxStepNum,m)

其中，zeros表示生成全0矩陣，maxStepNum為單個元胞包含單元最多數目；

步驟54，根據單元行列索引通過索引轉換得到單元在元胞-包含矩陣中的總索引，并通過索引賦值將對應編號單元填充進元胞-包含單元矩陣：

cellBallList(rowIndex+maxStepNum*(sortaIndex-1))＝sortI

其中，cellBallList表示元胞-包含單元矩陣；

步驟6，索引引用：

根據步驟54得到的元胞-包含單元矩陣，通過索引引用將步驟4得到的單元-鄰域元胞矩陣擴展為初步單元-鄰域單元矩陣：

nBallGrid0＝reshape(cellBallList(:,cIndex),27*maxStepNum,[])′

其中，cellBallList為元胞-包含單元矩陣，cIndex為單元-鄰域元胞矩陣，reshape與′為矩陣元素重排命令；

最終的初步單元-鄰域單元矩陣nBallGrid0為m行的矩陣，其每一行為相應單元鄰域元胞所包含單元的編號，不足部分值為0；按行降序排序：

nBallGrid0＝sort(nBallGrid0,2,′descend′)

此時值為0部分位于后列，去除全0列，記此時鄰域矩陣仍為nBallGrid0，剩余含0部分布爾索引：

zeroFilter＝nBallGrid0＝＝0

同時鄰域單元中包含自身單元編號，對應布爾索引：

selfFilter＝nBallGrid0＝＝(1:m)′

將自身單元及0單元用虛單元編號m+1填充,并重新排序:

nBallGrid0(zeroFilter|selfFilter)＝m+1

nBallGrid0＝sort(nBallGrid0,2)

去除自身單元引入的虛單元列：

nBallGrid＝nBallGrid0(:,1:(end-1))

nBallGrid即為最終單元-鄰域單元矩陣；

步驟7，距離篩選：基于虛單元填充方法實現坐標索引引用，進而通過距離篩選得到單元鄰居矩陣；

基于單元-鄰域單元矩陣，通過索引引用得到中心單元與鄰域單元的距離，沿p方向中心距離：

DP＝nP-IP

其中，nP為鄰域單元p方向坐標：

nP＝aP(nBallGrid)，aP包含虛單元坐標，

IP為中心單元p方向坐標：

IP＝P*ones(1,n0)

其中，ones為生成全一矩陣，n0為鄰居矩陣列數；

記x,y,z方向中心距離分別為DX,DY,DZ，距離判斷布爾矩陣：

DFilter＝max([DX,DY,DZ])<IR+nR+dSide

其中，IR為中心單元半徑，nR為鄰域單元半徑；

通過布爾矩陣篩選鄰居單元，并用虛單元編號m+1填充：

nBallGrid(～DFilter)＝m+1

其中，～表示布爾取反運算；

最后排序，去除多余虛單元列，即得單元-鄰居矩陣，記為nBall；

步驟8，基于鄰居矩陣的離散元求解計算；。