CN105760588B - 一种基于二层规则网格的sph流体表面重建方法 - Google Patents

Abstract

本发明提供一种基于二层规则网格的SPH流体表面重建方法，该方法基于CUDA平台并行实现，将整个模拟区域进行网格边长为2r的L1层粗粒度网格划分，待确定该层中流体表面附近区域的网格后，进一步细分成L2层精细网格；具体包括以下步骤：S1：确定模拟区域表面网格顶点；S2：计算标量值；S3：网格三角化。本发明的一种二层规则网格空间数据结构，提高了网格顶点在显存中存储的空间局部性，通过减少数据在显存和计算单元间的传输次数，大大提高了表面重建算法的整体运算效率,达到实时、高效渲染的目的。

Description

一种基于二层规则网格的SPH流体表面重建方法

技术领域

本发明属于流体现象模拟技术领域，具体涉及一种基于二层规则网格的SPH流体表面重建方法。

背景技术

光滑粒子流体动力学SPH(Smoothed Particle Hydrodynamics)模型被广泛应用到不同的流体现象模拟中，例如烟雾、泡沫、沙粒及基于大规模数据的复杂流体现象模拟等。它将连续的流体介质用离散的粒子表示，每个粒子承载各种物理量，包括质量，密度，速度等，通过跟踪每个粒子的运动轨迹，获取流体整体的运动行为。

在计算机图形学中，为了达到逼真的渲染效果，通常需要对基于粒子的流体模拟结果进行表面重建，即将离散的粒子进行表面提取，并用连续的网格进行表示。Lorensenand Cline提出的Marching Cubes(MC)方法是一种经典的等值面提取算法，它以整个模拟空间的三维标量域作为输入进行流体表面的提取，而三维标量域的计算及存储方式又对最终构网质量及效率产生很大影响。根据不同的网格划分方式通常分为规则网格和自适应网格两种构网方法。规则网格方法将整个模拟区域划分成规则网格，计算并存储每一个网格顶点的标量值，如Müller、Zhu and Bridson、Solenthaler等人的方法。该类方法优点是实现简单快捷，并可通过并行计算的方法提高效率，但由于最终提取表面相对于整个空间只是很薄的一层，该方法不可避免造成存储和计算上的冗余。自适应网格方法在流体表面附近区域用精细网格表示，而远离表面的区域用粗糙网格表示，常用方法有八叉树等，如Strain、Frisken等人的方法。该类方法可降低存储空间的占用，减少不必要网格顶点标量值的计算，但算法的并行实现较困难，同时八叉树等数据结构的构造、更新及查询算法复杂，对整体效率产生影响。

目前一种较为先进的方法是Solenthaler等人提出的基于规则网格划分的并行narrow-band方法。该方法将模拟区域划分成规则网格，相比于传统的计算并存储所有网格点的标量值，该方法只提取出流体表面附近区域的网格点进行计算，整个算法的框架基于并行方法实现。

普通规则网格法虽然实现简单，但在存储和运算效率上都存在缺陷，本发明正是基于这些缺陷进行改进，达到了提高存储和运算效率的目的；

另一种基于规则网格进行划分的是自适应多层网格法(例如Akinci等人的三层网格法)，该方法根据流体表面的曲率决定采用的网格级别，曲率大的采用更精细的网格，曲率小的采用较粗糙的网格，这在保证构网质量的基础上，提高了存储效率。但是，该方法由于存在不同网格层次间三角形的衔接问题，很难用并行方式实现，因此难以获取较高的运算效率。

相比于规则网格法，八叉树法也能达到自适应构建网格以降低存储占用的目的，但是八叉树数据结构本身的构造及结点查询复杂度高，并且也难以用并行方式实现，与流体表面重建的整体并行框架不能很好的结合，因此目前较少采用此方法。

发明内容

为降低流体表面构网的存储消耗、提高运算效率，本发明提供了一种基于二层规则网格的SPH流体表面重建方法，在保证构建高质量的流体表面网格基础上，进一步优化算法的存储和计算效率，实现针对大数据量流体粒子的表面高效准确渲染。

本发明的具体技术方案如下：

一种基于二层规则网格的SPH流体表面重建方法，该方法基于CUDA平台并行实现，将整个模拟区域进行网格边长为2r的L1层粗粒度网格划分，其中r为流体粒子间平衡距离；待确定该层中流体表面附近区域的网格后，进一步细分成L2层精细网格；具体包括以下步骤：具体包括以下步骤：

S1：确定模拟区域表面网格顶点；

S2：计算标量值；

S3：网格三角化。

其中，步骤S1的具体实施步骤如下：

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

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

其中d＝2r；待确定所有表面粒子占据的网格单元后，并行对每个网格搜寻相邻的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},Pi＝{P_j ⁱ|j＝0,1,…63}，CoarseSurfaceCells数组中位置i处对应元素值为C_i，则第i个线程按照如下规则将Vⁱ值写到数组位置Pⁱ中：

V_j ⁱ＝64·C_i+j (7)

P_j ⁱ＝64·i+j (8)

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

(3)哈希表构建：

构建哈希表只对L1层表面网格顶点进行哈希映射；数组CoarseSurfaceCells中每个元素作为键，对应的数组位置作为值构造键值对；键值对的哈希映射采用并行布谷鸟哈希(cuckoo hash)算法(该算法是Alcantara等于2009年在ACM Transactions on Graphics期刊上提出)，决定每个键值对的映射方式。

进一步，步骤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的标量值采用Solenthaler等人提出的特征值分析方法进行插值计算，表面隐函数定义如下：

其中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；用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算法(Lorensen和Cline于1987年在ACM SiggraphComputer Graphics杂志上提出)计算出要提取的三角形。

本发明的有益效果如下：

第一：基于CUDA平台并行实现流体粒子的表面重建，其中数据在显存和计算单元间的传输对算法整体效率有着重要的影响，本发明设计了一种二层规则网格空间数据结构，提高了网格顶点在显存中存储的空间局部性，通过减少数据在显存和计算单元间的传输次数，大大提高了表面重建算法的整体运算效率。具体为：二层规则网格空间数据结构将空间上连续的网格顶点在显存中尽可能连续存储，这样当搜索网格顶点的邻域粒子时，粒子在数组中也更大可能连续存储，这符合CUDA中显存的合并访问优化策略，可以以更少的读取次数获取所需要的数据，减少了数据传输耗费的时间。

第二、由于流体表面重建算法中所提取的表面只占整个模拟区域较小的比例，传统方法中分配与整个模拟区域相对应的存储空间会造成存储上的浪费，本研究采用哈希的方式只将与表面区域相关的信息进行映射，使得存储空间占用与模拟空间大小无关，极大的提高了存储效率。

第三、哈希算法选用并行布谷鸟算法实现，与本研究表面重建算法的并行框架有很好的的适应性，保证了哈希映射的运算效率。

附图说明

图1：基于二层规则网格结构的表面重建算法流程图；

图2：细分L1层网格单元二维示意图；

图3：网格顶点邻域粒子搜索示意图，其中中心顶点为待求解标量值的网格顶点，周围半径为2r的影响区域内所有粒子参与网格顶点标量值的插值求解；

图4：实施例中角落溃坝表面重建模拟初始状态示意图；

图5：实施例中角落溃坝表面重建模拟的中间状态示意图；

图6：实施例中角落溃坝表面重建模拟的最终状态示意图；

图7：实施例中索引排序法数据结构示意图。

具体实施方式

下面将结合本发明实施例，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

现有的流体表面构网方法以提高流体表面构网质量、降低存储消耗、提高运算效率为目标。但三者之间较难达到很好的平衡，比如提高网格精细度可以保留更多的表面细节，产生较高质量的流体表面，但相应增加了存储消耗和运算时间。Solenthaler等人提出的基于规则网格划分的并行narrow-band方法可以以较高的运算效率、较低的存储提取出高质量的流体表面。该方法通过建立整个模拟空间的三维标记数组标记出流体表面附近的区域，网格顶点标量值的计算和存储只与该区域的大小有关，而不受整个模拟空间的影响。虽然很大程度上提高了存储和运算效率，但三维标记数组本身仍然需要依据整个模拟空间规模而建，这在网格精度较高时仍会占用较大的存储空间。另外，算法通常未考虑网格顶点在显存中的存储顺序，而这对显存数据的读取速度有着至关重要的影响。

本发明在保证构建高质量的流体表面网格基础上，进一步优化算法的存储和计算效率，实现针对大数据量流体粒子的表面实时渲染。为了更好的对本发明的方法进行说明，本部分结全一个流体模拟中的经典场景角落溃坝进行表面重建的具体实施例进行说明。

实施例1

本案例对流体模拟中的经典场景角落溃坝进行表面重建。该场景模拟水池角落的巨大水块自由释放后的运动过程。本发明中算法基于任意粒子方法(如SPH)的模拟结果进行表面重建，以模拟结果的粒子位置、粒子速度作为算法输入。

场景参数配置：

模拟场景坐标方向参照OpenGL坐标系统。模拟区域范围为最小值MinBoundary(-250,0,-250)，最大值MaxBoundary(250,420,250)，在x和y方向的长度分别为L＝500，M＝420。水块由1.7M粒子进行初始化，水块范围为最小值MinParticles(-250,0,-250)，最大值MaxParticles(80,240,80)，粒子间距r＝2.212。本发明设计的二层规则网格空间数据结构中L1层网格单元边长为2r＝4.424，细分后L2层网格单元边长为r/2＝1.106。L1层网格在x方向的单元数目CellX＝L/2r≈114，在y方向的单元数目CellY＝M/2r≈95。

参照图1，本实施例具体实施步骤如下：

预处理步骤：首先采用Müller的光滑色域方法(smoothed color field method)标记出表面粒子作为预处理步骤。需分配与粒子数组对应的新数组，使用标志位bool类型进行标记。

第一步：确定表面网格顶点

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

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

C＝C3.x+C3.y·114+C3.z·114·95

参照图2，待确定所有表面粒子占据的网格单元后，并行对每个网格搜寻相邻的3³个L1层网格，作为该粒子所确定的L1层表面网格。标量值的计算在所确定的L1层表面网格区域进行。存储该L1层网格顶点的数组经过去重处理，得到最终确定的L1层表面网格顶点数组coarseSurfaceCells。

(2)细分L1层表面网格

所确定的L1层表面网格被进一步细分成64个小网格，作为L2层网格，每个网格单元边长为r/2＝1.106。分配新数组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ⁱ中：

V_j ⁱ＝64·C_i+j

P_j ⁱ＝64·i+j

(3)哈希表构建

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

第二步：标量值计算

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

将数组FineSurfaceVertices中每个表面网格顶点索引值除以64获取商P1和余数P2，再根据公式(4)和(5)分别反推出向量值CPos和IPos，带入公式(3)可求得该网格顶点的位置。

例如，数组FineSurfaceVertices中某个网格顶点索引值为1127，则可反推出P1＝17，P2＝39，进而算出CPos＝(17,0,0),IPos＝(3,1,2),FineVertexPos＝(71,1,2)。由此可计算出该顶点的坐标值：

PointCoordinates＝MinBoundary+r/2〃FineVertexPos＝(-171.474,1.106,-247.788)。

(2)搜寻邻域粒子

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

首先将整个模拟区域划分成规则网格，网格单元大小2r＝4.424。每个粒子对应一个网格，每个网格单元的索引值按照公式(1)(2)求解。粒子数组根据对应的网格索引进行排序，排序算法选用CUDA库提供的并行根排序算法，得到排序后粒子数组sortedParticlesArray。分配另一个数组GridCells，与划分的规则网格相对应，数组大小为所有网格单元数目，用于存储句柄对插入的粒子进行管理。每个粒子非空网格单元指向sortedParticlesArray数组中第一个插入的粒子。

用图7进一步说明索引排序算法过程：上层代表规则网格数组GridCells，网格内数字代表网格索引值，下层代表排序后粒子数组sortedParticlesArray。每个规则网格单元中若存在粒子，则指向数组sortedParticlesArray中的首个插入粒子。

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

查找到网格顶点邻域内所有粒子后，对于任意网格顶点位置x的标量值采用Solenthaler等人提出的特征值分析方法进行插值计算，求得的标量值存储到数组FineSurfaceVertices中的对应位置。

第三步：网格三角化

对于FineSurfaceVertices数组中每个表面网格顶点，需要查询所在网格其它顶点在数组中的位置，以便读取出相对应的标量值。

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

A＝(i，j，k)

neiIPos＝IPos+A

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

Pflag＝flagX+(flagy＜＜1)+(flagz＜＜2)

其中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

neiPIdx＝FineVertIdx+ConstantOffset

(b)NEI类型

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

neiCCellIdx＝P1+ConstantOffset

其中d＝2r＝4.424。用neiCCellIdx作为键值读出已经建立的哈希表中对应的值i。用neiIPos可以计算出该顶点在相邻L1层网格单元内的位置：

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

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

neiPIdx＝64·i+j

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

为了优化CUDA并行处理过程，程序执行时将提取出的L2层表面网格顶点划分成若干block，每个bolck对应512个线程，这样每个block最多管理8个L1层网格中的512个L2层网格顶点。由于从显卡的全局存储器读取数据比共享存储器慢很多，相比于为每个NEI类型顶点执行哈希值获取操作的方法，本实施案例提前将7个邻接L1层网格顶点索引值作为键获取对应的值，并将这些值存储到共享存储器中以提升读取效率。这样对于每个L1层网格中的64个L2层网格顶点，从全局存储器中读取哈希值操作只需执行7次，大大减少了不必要的显存数据传输次数。

实例对比结果：

为了突出本发明的优势，通过实例对本算法与Solenthaler等人的规则网格法(UG)和Akinci等人的Narrow-band(NB)方法进行对比，对比结果见表1，所有数据均为每帧的平均结果。其中t_es表示提取表面网格顶点的时间，t_sf为标量值计算时间，t_tri为三角化时间，t_tol为每帧表面重建总时间。MEM为算法占用显存平均大小。通过对比可以看出，本发明的算法相比于UG和NB分别获得2倍和3倍的加速比，而存储占用上获得将近3倍和5倍的提升，进一步验证了本发明算法的有效性。

方法名称	t<sub>es</sub>[msec]	t<sub>sf</sub>[msec]	t<sub>tri</sub>[msec]	t<sub>tol</sub>[sec]	MEM[GB]
						UG方法	__	4573	186.1	4.86	3.56
NB方法	222.3	2929.2	77.9	3.33	1.8
						本发明方法	38.9	1425.1	83.7	1.63	0.71

表1 表面重建效率对比

应当理解，虽然本说明书按照实施方式加以描述，但并非每个实施方式仅包含一个独立的技术方案，说明书的这种叙述方式仅仅是为清楚起见，本领域技术人员应当将说明书作为一个整体，各实施例中的技术方案也可以经适当组合，形成本领域技术人员可以理解的其他实施方式。

Claims (1)

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算法计算出要提取的三角形。

Images