[发明专利]一种基于二层规则网格的SPH流体表面重建方法

权利要求书

1.一种基于二层规则网格的SPH流体表面重建方法，其特征在于：

该方法基于CUDA平台并行实现，将整个模拟区域进行网格边长为2r的L1层粗粒度网格划分，其中r为流体粒子间平衡距离；待确定该层中流体表面附近区域的网格后，进一步细分成L2层精细网格；具体包括以下步骤：S1：确定模拟区域表面网格顶点，S2：计算标量值，S3：网格三角化；其中：

所述步骤S1的具体实施步骤如下：

(1)提取L1层表面网格顶点：

并行遍历已确定表面粒子，计算每个粒子对应的L1层网格索引序号；对于任何粒子位置x＝(x,y,z)，对应网格索引C计算公式如下：

其中d＝2r，参数x_min、y_min、z_min分别代表模拟区域范围最小值，L、M分别代表在x和y方向的长度；待确定所有表面粒子占据的网格单元后，并行对每个网格搜寻相邻的3³个L1层网格，作为该粒子所确定的L1层表面网格；标量值的计算在所确定的L1层表面网格区域进行；存储该L1层网格顶点的数组经过去重处理，得到最终确定的L1层表面网格顶点数组coarseSurfaceCells；

(2)细分L1层表面网格：

所确定的L1层表面网格被进一步细分成64个小网格，作为L2层网格，每个网格单元边长为r/2；

每个L2层网格单元在显存中通过该网格左下角的网格顶点索引表示；假设任意L2层网格顶点位置FineVertexPos表示如下：

FineVertexPos＝4·CPos+IPos (3)

其中CPos为该顶点所在的L1层网格的位置，IPos为该顶点位于所在L1层网格的内部位置，IPos每个分量取值范围为[0,4)，则该顶点的索引P按如下公式计算：

P2＝IPos.x+4·IPos.y+16·IPos.z (5)

P＝64·P1+P2 (6)

其中d＝2r；分配新数组FineSurfaceVertices，用于存储所有L2层网格顶点索引和对应的标量值，数组大小为数组coarseSurfaceCells大小的64倍；该存储过程通过并行方式处理，每个线程对应L1层网格中每个表面网格，并负责将该网格内的L2层网格顶点按照定义的索引计算方式存储到coarseSurfaceCells数组连续的位置；也就是说，假设Vⁱ＝{V_jⁱ|j＝0,1,…63},Pⁱ＝{P_jⁱ|j＝0,1,…63}，CoarseSurfaceCells数组中位置i处对应元素值为C_i，则第i个线程按照如下规则将Vⁱ值写到数组位置Pⁱ中：

由此L1层每个表面网格内空间连续的L2层网格顶点，在显存中也连续存储，但不同的L1层表面网格点在显存中却不一定连续存储；由于后续的网格三角化算法中需要确定L2层表面网格单元8个顶点在数组中的位置，对于8个顶点跨越不同L1层网格的情形，需要借助哈希表结构完成；

(3)哈希表构建：

构建哈希表只对L1层表面网格顶点进行哈希映射；数组CoarseSurfaceCells中每个元素作为键，对应的数组位置作为值构造键值对；键值对的哈希映射采用并行布谷鸟哈希算法，决定每个键值对的映射方式；

所述步骤S2的具体实施步骤如下：

标量值的计算在所提取的表面网格顶点上进行，基于数组中存储的表面网格顶点索引值计算网格顶点坐标，搜寻邻域内粒子以及进行标量值的插值计算，具体按以下步骤进行：

(1)计算L2层表面网格顶点坐标：

由于数组FineSurfaceVertices中每个表面网格顶点索引值由公式(6)求得，因此将该值除以64获取商P1和余数P2，再根据公式(4)和(5)分别反推出向量值CPos和IPos，带入公式(3)可求得该网格顶点的位置；根据预定义的模拟区域最小坐标值和网格单元大小r/2，可以求得该网格顶点的坐标值；

(2)搜寻邻域粒子：每个网格顶点周围影响半径为2r距离的区域作为邻接区域，该区域内所有粒子要对网格顶点的标量值进行插值计算；为了加速对粒子的邻接搜索，采用索引排序方法进行；

索引排序方法将整个模拟区域划分成规则网格，网格单元大小为2r；每个粒子对应一个网格，每个网格单元的索引值按照公式(1)、(2)求解；粒子数组根据对应的网格索引进行排序，得到排序后粒子数组sortedParticlesArray；分配另一个数组GridCells，与划分的规则网格相对应，数组大小为所有网格单元数目，用于存储句柄对插入的粒子进行管理；每个粒子非空网格单元指向sortedParticlesArray数组中第一个插入的粒子；该数据结构信息构造后，给定任意网格单元索引，可以查询出指向的排序后粒子数组中连续的粒子信息；

(3)计算网格顶点标量值：

查找到网格顶点邻域内所有粒子后，对于任意网格顶点位置x的标量值进行插值计算，表面隐函数定义如下：

其中r为粒子半径，R＝4r，W为密度核函数；因数f计算公式如下：

其中，EV_max为的最大特征值；求得的标量值存储到数组FineSurfaceVertices中的对应位置；

所述步骤S3的具体实施步骤如下：

本阶段采用MC方法并行对L2层每个表面网格进行三角划分，每个网格8个顶点的标量值作为MC方法的输入；对于FineSurfaceVertices数组中每个表面网格顶点，需要查询所在网格其它顶点在数组中的位置，以便读取出相对应的标量值；

对于数组FineSurfaceVertices中的任意表面网格顶点FineVertIdx，采用类似上述标量值计算中第一步计算L2层表面网格顶点坐标中的方法求得P1和IPos；该顶点所在网格的8个顶点位置neiIPos可以表示成：

A＝(i，j，k) (11)

neiIPos＝IPos+A (12)

其中矢量A的分量i，j，k取值0或1；如果求得的neiIPos的任意分量大于3，则表明该顶点属于相邻的L1层网格单元，标记为NEI类型；否则，该顶点属于同一网格单元；为了确定顶点的类型，设定标志位Pflag，按如下方式计算：

Pflag＝flagX+(flagY＜＜1)+(flagZ＜＜2) (13)

其中flagX，flagY，flagZ对应neiIPos的每个分量，且分量值大于3时赋值1，否则赋值0；Pflag取值范围为[0，8)，其中0表示顶点为INSIDE类型，其余值表示NEI类型，对应7个邻接L1层网格单元；

(a)INSIDE类型

该类型的顶点在数组FineSurfaceVertices中的数组位置neiPIdx按如下公式计算：

ConstantOffset＝A.x+4·A.y+16·A.z (14)

neiPIdx＝FineVertIdx+ConstantOffset (15)

(b)NEI类型该类型的顶点首先计算其所在L1层网格单元的索引neiCCellIdx：

neiCCellIdx＝P1+ConstantOffset (17)

其中d＝2r，L、M分别代表在x和y方向的长度；用neiCCellIdx作为键值读出已经建立的哈希表中对应的值i；用neiIPos可以计算出该顶点在相邻L1层网格单元内的位置neiIPos′：

neiIPos′＝(neiIPos.x％4，neiIPos.y％4，neiIpos.z％4) (18)

内部索引采用公式(5)求得，标记为j；最终该顶点的数组位置neiPIdx表示为：

neiPIdx＝64·i+j (19)

根据得到顶点在数组FineSurfaceVertices中的位置，读出对应的标量值；得到网格单元8个顶点的标量值后，通过MC算法计算出要提取的三角形。