本发明涉及网格处理技术领域，公开了一种分布式非结构网格跨处理器面对接方法及系统，该对接方法，采用两级索引结构并行化识别网格分区边界两侧的面对接关系；对任意两个对接面元素，依次对它们的质心坐标和归一化的网格点序列执行等价性判断。包括以下步骤：S1，并行导入分布式非结构网格基本几何数据；S2，跨处理器构造子区域之间的对偶通信列表；S3，并行构造对接边界处的曲面网格离散结构；S4，在各处理器上构造查询集合族的叉树结构；S5，跨处理器查询对接面元素之间的配对关系；S6，并行导出分布式非结构网格的面对接信息。本发明解决了现有技术存在的大规模非结构网格处理时处理效率低、数据处理能力差等问题。
1.一种分布式非结构网格跨处理器面对接方法，其特征在于，采用两级索引结构并行化识别网格分区边界两侧的面对接关系；
对任意两个对接面元素，依次对它们的质心坐标和归一化的网格点序列执行等价性判断；
S11，各处理器依次导入第0至M‑1级关联子网格的空间点信息；
S12，各处理器依次导入第0至M‑1级关联子网格的面元素信息；
S13，各处理器依次导入第0至M‑1级关联子网格的体单元信息；
S14，各处理器依次导入第0至M‑1级关联子网格上每个面元素两侧的体单元信息以及面元素的边界性质标识；当某个面元素位于子网格内部时，该面元素的两侧都有体单元存在；当某个面元素的边界性质标识显示该面元素是物理边界或区域分解边界时，且该面元素其中的一侧在当前处理器上不存在体单元，在该面元素的不存在体单元的一侧设置虚拟单元容器；
S21，计算全局一致的网格尺度公差epsilon：在各处理器上，对第0级子网格的面元素序列以及每个面元素的棱边进行二重遍历，设第i号面元素的第j号棱边的长度为length[i][j]，则当前计算区域上子网格的最小分辨率resolution=min{length[i][j],i=0,1,
2,…,j=0,1,2,…}，通过MPI全局规约操作得到resolution，然后利用epsilon=resolution*beta得到epsilon，再用epsilon来实现空间点的重合判断；其中，0.01≤beta≤0.1；
S22，构造第0级子网格坐标序列的三阶协方差矩阵G并计算协方差矩阵G的特征向量，三阶协方差矩阵G的计算公式为：G[i][j]=Covariance(g[i],g[j])；
其中，Covariance表示二元向量的协方差运算符，G[i][j]表示三阶协方差矩阵G的第i行第j列的元素；
在每个处理器上，设第0级子网格空间点坐标数据依次存储于三个K维向量容器中，Tg[0]=(x[0],x[1],…x[s],…,x[K‑1]) ,T
g[1]=(y[0],y[1],…x[s],…,y[K‑1]) ,T
g[2]=(z[0],z[1],…x[s],…,z[K‑1]) ，其中，s表示第0级子网格空间点的网格点编号，x[s],y[s],z[s]为第s个网格点的三个空间坐标分量，K表示网格点序列的长度，T表示向量的转置运算；通过数值方法计算协方差矩阵G的全部单位特征向量e[0]、e[1]、e[2]，并使单位特征向量e[0]、e[1]、e[2]满足右手系规则；
S23，在各处理器上构造第0级子网格的最小边界包围盒；
根据步骤S22的结果，在三维欧氏空间中选取一个长方体，此长方体同时满足以下条件：①三组平行的表面分别以单位特征向量e[0]、e[1]、e[2]为法向量；②第0级子网格完全包含在长方体区域内部；③子网格到长方体六个表面的最小距离皆等于全局一致的网格尺度公差epsilon；以此长方体作为当前处理器上第0级子网格的最小边界包围盒，结合MPI广播通信，使每个最小边界包围盒的几何参数在全部处理器上实现共享存储；
S24，在各处理器上并行判断最小边界包围盒之间的相交关系，建立子区域之间的对偶通信二维表。
2.根据权利要求1所述的一种分布式非结构网格跨处理器面对接方法，其特征在于，步骤S3包括以下步骤：S31，在各处理器上，根据边界性质标识依次提取第m级子网格上对接面元素的局部编号，分别存储在C++标准模板库的向量容器中并关联到相应子网格上；
S32，在各处理器上，依次提取第m级子网格上对接面元素包含的网格点总数，分别存储在C++标准模板库的向量容器中并关联到相应子网格上；
S33，在各处理器上，根据S11的空间点信息和S12的面元素信息建立面元素与网格点之间的关联关系，依次提取第m级子网格上对接面元素具备的网格点坐标数据，分别存储在C++标准模板库的向量容器中并关联到相应子网格上。
3.根据权利要求2所述的一种分布式非结构网格跨处理器面对接方法，其特征在于，步骤S4包括以下步骤：S41，在各处理器上，对第m级子网格的对接面容器，以S21计算出的全局一致的网格尺度公差epsilon作为对接面网格点序列的交替方向数字树latticeBinaryTree和对接面质心序列的交替方向数字树centroidBinaryTree的公差参数，完成latticeBinaryTree和centroidBinaryTree的根结点的初始化；
S42，在各处理器上，对第m级子网格的对接面容器，依次遍历每个对接面上的网格点序列，生成latticeBinaryTree的各个子结点，同时完成网格点的归一化；
S43，在各处理器上，对第m级子网格的对接面容器，依次遍历每个对接面上的质心序列，生成centroidBinaryTree的各个子结点，完成质心等价面元素的初筛分组。
4.根据权利要求3所述的一种分布式非结构网格跨处理器面对接方法，其特征在于，步骤S43中，完成质心等价面元素的初筛分组的过程中，如果在centroidBinaryTree上不存在包含当前对接面元素的质心等价分组，则在相应位置上增加新的子结点，并与新建的等价对接面地址空集合链接起来，然后，把当前对接面元素的首地址插入该集合中；如果在centroidBinaryTree存在一个子结点，与之链接的对接面地址集合包含与当前对接面质心等价的元素，则无需增加新的子结点，直接将前对接面的地址插入到该集合中。
5.根据权利要求4所述的一种分布式非结构网格跨处理器面对接方法，其特征在于，步骤S5包括以下步骤：S51，在第n号处理器上，对第m级子网格的对接面元素序列，将每个对接面的编号、网格点总数以及空间坐标数据加以打包，发送到对偶处理器上进行信息共享；
S52，在每个对偶处理器上，依次解析打包数据中的对接面信息，计算对接面元素的质心坐标，并结合同级子网格对应的交替方向数字树centroidBinaryTree加以查询比较，在latticeBinaryTree上找到全部的关联网格点；
S53，在对偶处理器上，结合同级子网格相关的交替方向数字树latticeBinaryTree，对当前解析出的对接面元素网格点序列进行编号并排序，然后进入质心等价分组中执行精确查询与比较；如果质心等价分组中存在各网格点编号完全一致的元素，说明当前面对接元素配对成功，记录相应的对偶元素编号和当前对偶处理器的编号，否则说明配对失败，直接删除当前来自其它处理器的对接面元素的几何信息；
S54，在第n号处理器的一个或多个对偶处理器上，记录第m级子网格上配对成功的几何信息并打包，结合MPI通信函数接口，逆向返回到第n号处理器上加以解析，将面对接关系存储在S14构造的虚拟单元容器中。
6.一种分布式非结构网格跨处理器面对接系统，其特征在于，基于权利要求1至5任一项所述的一种分布式非结构网格跨处理器面对接方法，包括依次电相连的并行导入模块、对偶通信构建模块、离散结构构建模块、叉树结构构建模块、配对关系查询模块、并行导出模块，并行导入模块，用以并行导入分布式非结构网格基本几何数据；
对偶通信构建模块，用以跨处理器构造子区域之间的对偶通信列表；
离散结构构建模块，用以并行构造对接边界处的曲面网格离散结构；
叉树结构构建模块，在各处理器上构造查询集合族的叉树结构；
配对关系查询模块，跨处理器查询对接面元素之间的配对关系；
并行导出模块，并行导出分布式非结构网格的面对接信息。
一种分布式非结构网格跨处理器面对接方法及系统
技术领域
[0001] 本发明涉及网格处理技术领域，具体是一种分布式非结构网格跨处理器面对接方法及系统。
背景技术
[0002] 非结构网格技术是当今计算流体力学软件对复杂构型内外流场进行数值模拟的重要依托，它不仅能够在非规则区域上进行灵活的单元离散，大幅降低网格生成的时间成本，而且可以通过图切割算法实现快速的区域分解，达到高度的负载平衡以适应现代高性能计算的需求。目前广泛采用的METIS与PARAMETIS等计算几何工具库已成功实现图分解算法的封装，用户仅需调用标准的函数接口即可迅速确定任意计算单元所属的处理器编号。在此基础上，MGRIDGEN等计算软件能够进一步执行体单元聚合操作，构造出几何多重网格的空间数据结构，以满足求解器在物理场加速收敛方面的高级需求。
[0003] 然而，仅依靠上述计算几何工具库并不能得到网格在分区边界两侧的面对接关系，开发人员仍需要编写大量的代码完成复杂的查询运算，以此建立大规模并行计算所需的通信映射。伴随着网格规模的提升，每个子区域及其相邻区域上的对接面元素总量均同步增长，从而使关联映射的总时间依非线性规律延长。同时，现有的关联映射算法几乎都是串行实现的，不能深度挖掘高性能集群的计算潜力。以国产CFD软件PHengLEI和FlowStar为例，当非结构网格规模达到十亿量级时，其处理相邻子区域面对接关系的内核算法执行时间高达数小时甚至数天。这些因素，成为大规模非结构网格在高性能计算方面的严重障碍。
[0004] 现有技术存在如下问题：分布式非结构网格跨处理器通信关系难以并行化建立，从而使几何前处理过程产生大量的时间开销，网格处理能力难以适应大规模并行计算的需求。
发明内容
[0005] 为克服现有技术的不足，本发明提供了一种分布式非结构网格跨处理器面对接方法及系统，解决现有技术存在的大规模非结构网格处理时处理效率低、数据处理能力差等问题。
[0006] 本发明解决上述问题所采用的技术方案是：
[0007] 一种分布式非结构网格跨处理器面对接方法，采用两级索引结构并行化识别网格分区边界两侧的面对接关系。
[0008] 作为一种优选的技术方案，对任意两个对接面元素，依次对它们的质心坐标和归一化的网格点序列执行等价性判断。
[0009] 作为一种优选的技术方案，包括以下步骤：
[0010] S1，并行导入分布式非结构网格基本几何数据；
[0011] S2，跨处理器构造子区域之间的对偶通信列表；
[0012] S3，并行构造对接边界处的曲面网格离散结构；
[0013] S4，在各处理器上构造查询集合族的叉树结构；
[0014] S5，跨处理器查询对接面元素之间的配对关系；
[0015] S6，并行导出分布式非结构网格的面对接信息。
[0016] 作为一种优选的技术方案，步骤S1包括以下步骤：
[0017] S11，各处理器依次导入第0至M‑1级关联子网格的空间点信息；
[0018] S12，各处理器依次导入第0至M‑1级关联子网格的面元素信息；
[0019] S13，各处理器依次导入第0至M‑1级关联子网格的体单元信息；
[0020] S14，各处理器依次导入第0至M‑1级关联子网格上每个面元素两侧的体单元信息以及面元素的边界性质标识；当某个面元素位于子网格内部时，该面元素的两侧都有体单元存在；当某个面元素的边界性质标识显示该面元素是物理边界或区域分解边界时，且该面元素其中的一侧在当前处理器上不存在体单元，在该面元素的不存在体单元的一侧设置虚拟单元容器。
[0021] 作为一种优选的技术方案，步骤S2包括以下步骤：
[0022] S21，计算全局一致的网格尺度公差epsilon：在各处理器上，对第0级子网格的面元素序列以及每个面元素的棱边进行二重遍历，设第i号面元素的第j号棱边的长度为length[i][j]，则当前计算区域上子网格的最小分辨率resolution=min{length[i][j],i=0,1,2,…,j=0,1,2,…}，通过MPI全局规约操作得到resolution，然后利用epsilon=resolution*beta得到epsilon，再用epsilon来实现空间点的重合判断；其中，0.01≤beta≤0.1；
[0023] S22，构造第0级子网格坐标序列的三阶协方差矩阵G并计算协方差矩阵G的特征向量，三阶协方差矩阵G的计算公式为：
[0024] G[i][j]=Covariance(g[i],g[j])；
[0025] 其中，Covariance表示二元向量的协方差运算符，G[i][j]表示三阶协方差矩阵G的第i行第j列的元素；
[0026] 在每个处理器上，设第0级子网格空间点坐标数据依次存储于三个K维向量容器中，
[0027] g[0]=(x[0],x[1],…x[s],…,x[K‑1])T,
[0028] g[1]=(y[0],y[1],…x[s],…,y[K‑1])T,
[0029] g[2]=(z[0],z[1],…x[s],…,z[K‑1])T，
[0030] 其中，s表示第0级子网格空间点的网格点编号，x[s],y[s],z[s]为第s个网格点的三个空间坐标分量，K表示网格点序列的长度，T表示向量的转置运算；通过数值方法计算协方差矩阵G的全部单位特征向量e[0]、e[1]、e[2]，并使单位特征向量e[0]、e[1]、e[2]满足右手系规则；
[0031] S23，在各处理器上构造第0级子网格的最小边界包围盒；
[0032] 根据步骤S22的结果，在三维欧氏空间中选取一个长方体，此长方体同时满足以下条件：①三组平行的表面分别以单位特征向量e[0]、e[1]、e[2]为法向量；②第0级子网格完全包含在长方体区域内部；③子网格到长方体六个表面的最小距离皆等于全局一致的网格尺度公差epsilon；以此长方体作为当前处理器上第0级子网格的最小边界包围盒，结合MPI广播通信，使每个最小边界包围盒的几何参数在全部处理器上实现共享存储；
[0033] S24，在各处理器上并行判断最小边界包围盒之间的相交关系，建立子区域之间的对偶通信二维表。
[0034] 作为一种优选的技术方案，步骤S3包括以下步骤：
[0035] S31，在各处理器上，根据边界性质标识依次提取第m级子网格上对接面元素的局部编号，分别存储在C++标准模板库的向量容器中并关联到相应子网格上；
[0036] S32，在各处理器上，依次提取第m级子网格上对接面元素包含的网格点总数，分别存储在C++标准模板库的向量容器中并关联到相应子网格上；
[0037] S33，在各处理器上，根据S11的空间点信息和S12的面元素信息建立面元素与网格点之间的关联关系，依次提取第m级子网格上对接面元素具备的网格点坐标数据，分别存储在C++标准模板库的向量容器中并关联到相应子网格上。
[0038] 作为一种优选的技术方案，步骤S4包括以下步骤：
[0039] S41，在各处理器上，对第m级子网格的对接面容器，以S21计算出的全局一致的网格尺度公差epsilon作为对接面网格点序列的交替方向数字树latticeBinaryTree和对接面质心序列的交替方向数字树centroidBinaryTree的公差参数，完成latticeBinaryTree和centroidBinaryTree的根结点的初始化；
[0040] S42，在各处理器上，对第m级子网格的对接面容器，依次遍历每个对接面上的网格点序列，生成latticeBinaryTree的各个子结点，同时完成网格点的归一化；
[0041] S43，在各处理器上，对第m级子网格的对接面容器，依次遍历每个对接面上的质心序列，生成centroidBinaryTree的各个子结点，完成质心等价面元素的初筛分组。
[0042] 作为一种优选的技术方案，步骤S43中，完成质心等价面元素的初筛分组的过程中，如果在centroidBinaryTree上不存在包含当前对接面元素的质心等价分组，则在相应位置上增加新的子结点，并与新建的等价对接面地址空集合链接起来，然后，把当前对接面元素的首地址插入该集合中；如果在centroidBinaryTree存在一个子结点，与之链接的对接面地址集合包含与当前对接面质心等价的元素，则无需增加新的子结点，直接将前对接面的地址插入到该集合中。
[0043] 作为一种优选的技术方案，步骤S5包括以下步骤：
[0044] S51，在第n号处理器上，对第m级子网格的对接面元素序列，将每个对接面的编号、网格点总数以及空间坐标数据加以打包，发送到对偶处理器上进行信息共享；
[0045] S52，在每个对偶处理器上，依次解析打包数据中的对接面信息，计算对接面元素的质心坐标，并结合同级子网格对应的交替方向数字树centroidBinaryTree加以查询比较，在latticeBinaryTree上找到全部的关联网格点；
[0046] S53，在对偶处理器上，结合同级子网格相关的交替方向数字树latticeBinaryTree，对当前解析出的对接面元素网格点序列进行编号并排序，然后进入质心等价分组中执行精确查询与比较；如果质心等价分组中存在各网格点编号完全一致的元素，说明当前面对接元素配对成功，记录相应的对偶元素编号和当前对偶处理器的编号，否则说明配对失败，直接删除当前来自其它处理器的对接面元素的几何信息；
[0047] S54，在第n号处理器的一个或多个对偶处理器上，记录第m级子网格上配对成功的几何信息并打包，结合MPI通信函数接口，逆向返回到第n号处理器上加以解析，将面对接关系存储在S14构造的虚拟单元容器中。
[0048] 一种分布式非结构网格跨处理器面对接系统，基于所述的一种分布式非结构网格跨处理器面对接方法，包括依次电相连的并行导入模块、对偶通信构建模块、离散结构构建模块、叉树结构构建模块、配对关系查询模块、并行导出模块，
[0049] 并行导入模块，用以并行导入分布式非结构网格基本几何数据；
[0050] 对偶通信构建模块，用以跨处理器构造子区域之间的对偶通信列表；
[0051] 离散结构构建模块，用以并行构造对接边界处的曲面网格离散结构；
[0052] 叉树结构构建模块，在各处理器上构造查询集合族的叉树结构；
[0053] 配对关系查询模块，跨处理器查询对接面元素之间的配对关系；
[0054] 并行导出模块，并行导出分布式非结构网格的面对接信息。
[0055] 本发明相比于现有技术，具有以下有益效果：
[0056] （1）本发明在多处理器上构造出以嵌套叉树结构形成的集合族，通过二级查询方式深度挖掘高性能计算平台的数据处理能力，快速建立用于大型非结构网格跨处理器的面对接映射关系；
[0057] （2）本发明将大规模非结构网格的几何数据分布存储于多个处理器上（现有技术多为串行算法），有效缓解了几何前处理算法执行过程中面临的内存压力；
[0058] （3）本发明在以嵌套叉树结构形成的集合族上进行多级筛查，使目标对象的待检索范围迅速缩小（现有技术中面对接关系的判断依据单一），从而提升了不同子区域间面对接关系的识别效率。
附图说明
[0059] 图1是本发明所述的一种分布式非结构网格跨处理器面对接方法的步骤示意图。
[0060] 图2是本发明的分布式多重网格结构图。
[0061] 图3是本发明的网格面对接关系示意图。
[0062] 图4是本发明的子区域之间的对偶通信关系图。
[0063] 图5是本发明项目在每个处理器上存储的离散对接面数据结构图。
[0064] 图6是本发明第m级子网格对接面形成的集合族叉树结构。
[0065] 图7是本发明对接面元素跨处理器配对数据流示意图。
[0066] 图8为本发明实施中步骤S52、S53、S54的流程图。
具体实施方式
[0067] 下面结合实施例及附图，对本发明作进一步的详细说明，但本发明的实施方式不限于此。
[0068] 实施例
[0069] 如图1至图8所示，针对上述问题，本发明提供一种分布式内存环境下的非结构网格面对接技术，构造并行化的分析模型快速识别网格分区边界两侧的面对接关系。该模型的核心部分由“初筛+确认”两级索引结构组成。对任意两个对接面元素，分析模型依次对它们的质心坐标和归一化的网格点序列执行等价性判断。传统查询算法可以归纳为单处理器上两个集合的求交运算，相比之下，本发明在多处理器上构造出以嵌套叉树结构形成的集合族，通过二级查询方式深度挖掘高性能计算平台的数据处理能力，快速建立用于大型非结构网格跨处理器的面对接映射关系。
[0070] 为实现上述目的，根据附图1所示，本发明项目采用的技术解决方案为：一种分布式非结构网格跨处理器面对接技术。它的实施包含如下6个基本步骤（其中，S3至S5是一个循环过程，它们跨处理器识别出多重网格模型中各级子网格之间的面对接关系）：
[0071] S1：并行导入分布式非结构网格基本几何数据。当MPI环境初始化后，根据物理场模拟器中多重网格计算模型的需要，各处理器同步导入关联子区域上各级子网格的拓扑结构，并使空间体单元的分组编号与所在的处理器编号保持一致。本发明项目中的最稠密级子网格指多重网格的第0级，使用体单元聚合算法递推得到第1至M‑1级粗网格，其中M为多重网格分级数。附图2展示了本发明适用的分布式多重网格结构，在每个处理器上，箭头指向子网格逐级粗化的方向；虚线框内是全体处理器上同级子网格的分布式存储载体，本发明在每个虚线框内分别建立各级子网格跨处理器的面对接关系。
[0072] S2：跨处理器构造子区域之间的对偶通信列表。在整个求解区域上，每个子网格一般与多个邻居网格具有面对接关系。S1完成后，每个处理器上的独立子网格并不能自动感知各邻居子网格所在处理器的编号。因为每个处理器上的第0级子网格能够使子区域的边界形状保持不变，所以本发明项目为其构造最小边界包围盒，通过包围盒的相对位置识别各处理器之间的基本通信关系。
[0073] S3：并行构造对接边界处的曲面网格离散结构。结合附图2，本发明项目需在分区边界两侧处理同级子网格的面对接关系，一般情况下，在每个处理器上，对接面元素承载的几何数据并非分布在连续内存上。为了便于待搜索几何信息的跨处理器传输以及查询配对结果的返回，本发明项目使用独立的动态向量容器依次复制并存储对接面元素的几何信息。附图5是本发明项目在单个处理器上存储的离散对接面数据结构图，在每一级网格上，对接面元素从0开始编号，并且按照它在原始网格上的遍历顺序逐次增加。
[0074] S4：在各处理器上构造查询集合族的叉树结构。显然，空间两个直边面元素Face_A与Face_B重合的充要条件为：①Face_A与Face_B具备的网格点总数相等且关联网格点的空间位置相同。根据网格点序列的代数平均运算可得对接面元素的质心坐标，由此得出对接面元素Face_A与Face_B重合的一个必要条件：②Face_A与Face_B质心的空间位置相同。对于非结构网格面元素上的网格点而言，除三角形情形外，一般不能严格分布在同一平面内，而且当体单元聚合后，其表面多边形的扭曲更加明显，所以②是一个非充分条件。本发明项目以条件②为质心初筛判据对子网格上的对接面元素地址进行整体分组，并以叉树结构对各组地址进行排序存储。对每个面元素包含的空间网格点进行归一化后，结合C++标准模板库，使用基于整型数据的向量（vector）容器对其编号进行排序存储，然后在新的层次上使用集合（set）容器对这些质心等价的向量容器首地址进行排序存储。这样，第m级子网格对接面形成了附图6所示的集合族，它在计算机内部以嵌套形式的叉树结构加以存储，其中用实线圆圈表示的每个结点都对应于当前子网格上一个质心等价的对接面首地址集合。
[0075] S5：跨处理器查询对接面元素之间的配对关系。附图7是本发明对接面元素跨处理器配对过程中的数据流示意图，其中带双向箭头的虚线为查询操作前后信息的传输方向。在第n号处理器上，对第m级子网格，将对接面序列的几何信息加以顺序打包，通过MPI函数接口传输到对偶处理器上加以共享。在每个对偶处理器上，对上述打包的几何信息加以解析，复原它们的数据结构，并依托同级子网格对接面元素形成的嵌套叉树结构进行快速查询操作与配对判断。
[0076] S6：并行导出分布式非结构网格的面对接信息。
[0077] 进一步地，S1的具体操作步骤可以分解为：
[0078] S11：各处理器依次导入第0至M‑1级关联子网格的空间点信息，它们包含网格点序列的坐标数据以及该序列的长度。
[0079] S12：各处理器依次导入第0至M‑1级关联子网格的面元素信息，它们包含每个面元素上网格点序列在当前处理器上的局部编号以及该序列的长度。
[0080] S13：各处理器依次导入第0至M‑1级关联子网格的体单元信息，它们包含每个体单元上面元素序列在当前处理器上的局部编号以及该序列的长度。
[0081] S14：各处理器依次导入第0至M‑1级关联子网格上每个面元素两侧的体单元信息以及面元素的边界性质标识。当面元素位于子网格内部时，它的两侧都有体单元存在；当面元素的边界性质标识显示它是物理边界或区域分解边界时，其中一侧在当前处理器上不存在体单元，本发明项目在此设置虚拟单元，在后续步骤中，它们用来存储跨处理器的面对接信息。
[0082] 进一步地，S2的具体操作步骤可以分解为：
[0083] S21：计算全局一致的网格尺度公差epsilon。在各处理器上，对第0级子网格的面元素序列以及每个面元素的棱边进行二重遍历，设第i号面元素的第j号棱边的长度为length[i][j]，那么当前计算区域上子网格的最小分辨率resolution=min{length[i][j], i=0,1,2,…, j=0,1,2,…}，通过MPI全局规约操作即得到resolution的最小值，这里不妨仍以resolution代替，本发明项目规定全局一致的网格尺度公差epsilon=resolution*0.1，即比全局网格最小分辨率低一个数量级。epsilon为本发明项目关键的尺度参数，用来实现空间点的重合判断。
[0084] S22：构造第0级子网格坐标序列的协方差矩阵并计算其特征向量。在每个处理器上，设第0级子网格空间点坐标数据依次存储于三个K维向量容器中，即:
[0085] g[0]=(x[0],x[1],…x[s],…,x[K‑1])T,
[0086] g[1]=(y[0],y[1],…x[s],…,y[K‑1])T,
[0087] g[2]=(z[0],z[1],…x[s],…,z[K‑1])T，
[0088] 其中，x[s],y[s],z[s]为第s个网格点的三个空间坐标分量，K是网格点序列的长度。三阶协方差矩阵G的第i行第j列分量G[i][j]（i和j从0开始编号）的表达式为:
[0089] G[i][j]=Covariance(g[i],g[j]),
[0090] 这里Covariance为二元向量的协方差运算符。通过数值方法计算协方差矩阵G的全部单位特征向量e[0]，e[1]，e[2]，并使之满足右手系规则。
[0091] S23：在各处理器上构造第0级子网格的最小边界包围盒。根据S22的结果，在三维Euclid空间中选取一个长方体，它满足条件：①三组平行的表面分别以单位特征向量e[0]，e[1]，e[2]为法向量；②第0级子网格完全包含在长方体区域内部；③子网格到长方体六个表面的最小距离皆等于全局一致的网格尺度公差epsilon。这个长方体是存在且唯一的，以此作为当前处理器上第0级子网格的最小边界包围盒。第③个条件充分保证了任意子网格点是包围盒的内点，避免了点在边界处的复杂分析。结合MPI广播通信，使每个最小边界包围盒的几何参数在全部处理器上实现共享存储。附图3是本发明的网格面对接关系示意图，由于第0级子网格的最小边界包围盒是其外形的一个几何逼近，所以子网格之间的相邻关系可以通过最小边界包围盒的相交判断加以初步识别。
[0092] S24：在各处理器上并行判断最小边界包围盒之间的相交关系，建立子区域之间的对偶通信二维表。附图4是本发明的子区域之间的对偶通信关系图，其中虚线框内是对偶处理器的编号。以第1号处理器为例，结合附图3可知，当前子网格仅与第0、2、4、5、6号处理器上的子网格存在可能的对接关系。需要说明的是，本发明项目中，包围盒相交是子网格面对接的一个必要非充分条件。尽管如此，这个筛查条件有效避免了当前子网格与其它所有同级子网格之间不加区别的相邻判断，从而使跨处理器的通信成本得以显著降低。
[0093] 进一步地，S3的具体操作步骤可以分解为：
[0094] S31：在各处理器上，根据边界性质标识依次提取第m级子网格上对接面元素的局部编号，分别存储在C++标准模板库的向量（vector）容器中并关联到相应子网格上。
[0095] S32：在各处理器上，依次提取第m级子网格上对接面元素包含的网格点总数，分别存储在C++标准模板库的向量（vector）容器中并关联到相应子网格上。
[0096] S33：在各处理器上，根据S11和S12建立的面元素与网格点之间的映射关系，依次提取第m级子网格上对接面元素具备的网格点坐标数据，分别存储在C++标准模板库的向量（vector）容器中并关联到相应子网格上。
[0097] 进一步地，S4的具体操作步骤可以分解为：
[0098] S41：在各处理器上，对第m级子网格的对接面容器，以S21计算出的全局一致的网格尺度公差epsilon作为对接面网格点序列之交替方向数字树latticeBinaryTree和对接面质心序列之交替方向数字树centroidBinaryTree的公差参数，完成它们的根结点的初始化。
[0099] S42：在各处理器上，对第m级子网格的对接面容器，依次遍历每个对接面上的网格点序列，生成latticeBinaryTree的各个子结点，同时完成网格点的归一化，即空间坐标相同的网格点使用相同的几何对象加以记录。任意给定一个对接面元素，不妨设它的网格点总数为S，归一化之后的网格点编号序列为{p[0],p[1],…,p[S‑1]}，本发明项目首先使用C++标准模板库中的整型向量（vector）容器series记录这个编号序列，然后，结合泛型算法，按照编号值的大小对其排序，结果仍存储在series容器中。
[0100] S43：在各处理器上，对第m级子网格的对接面容器，依次遍历每个对接面上的质心序列，生成centroidBinaryTree的各个子结点，完成质心等价面元素的初筛分组。根据递归或非递归的查询结果将上述过程划分为两种情形：①如果在centroidBinaryTree上不存在包含当前对接面元素的质心等价分组，那么在相应位置上增加新的子结点，并与新建的等价对接面地址空集合链接起来，然后，把当前对接面元素的首地址插入该集合中。②如果在centroidBinaryTree存在一个子结点，与之链接的对接面地址集合包含与当前对接面质心等价的元素，那么无需增加新的子结点，直接将前对接面的地址插入到该集合中。需要说明的是，C++标准模板库中的集合（set）容器本质上也是一个红黑树排序结构，对S42定义的网格点编号序列，本发明项目使用C++标准模板库中的整型向量（vector）容器默认的排序规则构造函数对象，用于分组集合内部元素（即上述若干个vector对象）的排序运算。
[0101] 进一步地，S5的具体操作步骤可以分解为：
[0102] S51：在第n号处理器上，对第m级子网格的对接面元素序列，将每个对接面的编号、网格点总数以及空间坐标数据加以打包，发送到对偶处理器上进行信息共享。
[0103] S52：在每个对偶处理器上，依次解析打包数据中的对接面信息，计算其质心坐标，并结合同级子网格对应的交替方向数字树centroidBinaryTree（图8中简写为CBT）加以查询与比较。对于来自于其它处理器上对接面元素，上述查询过程划分为以下情形（如图8所示）：①判断GP是否存在：如果在centroidBinaryTree上不存在包含当前对接面元素的质心等价分组，那么在当前处理器上直接删除该对接面元素的几何信息（即图8中释放SF[s]的几何数据）。②然后，判断是否全部找到：如果在当前对接面元素上存在一点p，其坐标在遍历同级子网格对应的交替方向数字树latticeBinaryTree后仍无法查到重合点，那么在当前处理器上直接删除该对接面元素的几何信息（即图8中释放SF[s]的几何数据）；如果在centroidBinaryTree上存在当前对接面元素的质心等价分组，那么执行S53（本步骤中完成所有对接面元素的查询配对操作，找到所有的关联网格点）。
[0104] S53：参照S42方法，在对偶处理器上，结合同级子网格相关的交替方向数字树latticeBinaryTree，对当前解析出的对接面元素网格点序列进行编号并排序，然后进入质心等价分组（即若干个整型向量（vector）复合而成的集合（set）容器）中执行精确查询与比较。上述过程由C++标准模板库中集合（set）内置的红黑树快速查询算法加以实现。如果质心等价分组中存在各网格点编号完全一致的元素，说明当前面对接元素配对成功，记录相应的对偶元素编号和当前对偶处理器的编号，否则说明配对失败，直接删除当前来自其它处理器的对接面元素的几何信息。
[0105] S54：在第n号处理器的一个或多个对偶处理器上，记录第m级子网格上配对成功的几何信息并打包，结合MPI通信函数接口，逆向返回到第n号处理器上加以解析，将面对接关系存储在S14构造的虚拟单元容器中。
[0106] 如上所述，可较好地实现本发明。
[0107] 本说明书中所有实施例公开的所有特征，或隐含公开的所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以以任何方式组合和/或扩展、替换。
[0108] 以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，依据本发明的技术实质，在本发明的精神和原则之内，对以上实施例所作的任何简单的修改、等同替换与改进等，均仍属于本发明技术方案的保护范围之内。
