一种基于辅助树的多层快速多极子并行网格细剖方法

摘要

本发明公开了一种基于辅助树的多层快速多极子并行网格细剖方法，本发明基于粗网格三角形单元中心构建辅助树，对此辅助树采取与计算所用多极子树结构相同的离散模式、进程数进行分层离散，并根据按盒子并行首层的离散模式，经遍历到最细层建立三角形单元与计算所用进程间的分布映射，对三角形单元、边、顶点编号重排形成分区连续的粗网格文件；其次，并行读取粗网格文件，各进程上对粗网格进行均匀一致性细剖，并对新生成的边、点、三角形重编号形成完整细网格信息，本发明方法可最大限度的保持数据的本地性，缩减通信数据量，因此具有很高的并行效率。

权利要求

1.一种基于辅助树的多层快速多极子并行网格细剖方法，其特征在于，包括以下步骤：步骤一、将复杂结构目标剖分为三角形网格，根据三角形网格的顶点编号和顶点坐标计算三角形中心坐标，基于三角形中心，按照多层快速多极子盒子划分原则，生成网格的辅助树结构；

步骤二、根据预设的并行方案和MPI进程数，对网格进行负载分配，确定归属于各进程上的三角形网格编号；

步骤三、对分配后的三角形网格编号、顶点编号和棱边编号进行重排序，并输出对应的存储文件；

对所述三角形网格编号、顶点编号和棱边编号进行重排序具体为：三角形网格编号：首先对离散到进程上不连续的三角形单元编号按照原编号大小顺序压缩重排为自1开始的连续编号，之后按照进程增序连续编号；

顶点编号：先将所有进程内的顶点编号，然后将公共顶点按顺序重新编号；

棱边编号：对进程内的棱边按原编号压缩重排后按照进程增序编号，然后对所有公共边压缩重排后编号，公共边编号的起始位置在所有进程内边编号之后；

步骤四、并行读取上述步骤三得到的存储文件，对三角形网格进行预设次数的细剖，并对新生成的棱边、顶点和三角形网格重排序，形成完整的细网格信息；

对所述新生成的棱边、顶点和三角形网格重排序，具体为：三角形网格编号和顶点编号采用所有进程的全局编号；

顶点编号的原有顶点编号值不变，新产生的顶点按照从属边编号压缩后局部顺序重新编号。

2.如权利要求1所述的一种基于辅助树的多层快速多极子并行网格细剖方法，其特征在于，步骤四中，细剖的剖分尺度m与次数按如下方式确定：m≤m0×2n

其中，n是预设的细剖次数，m0是最终需要达到的平均剖分边长，设置为十分之一波长。

说明书

一种基于辅助树的多层快速多极子并行网格细剖方法

技术领域

[0001] 本发明属于计算电磁计算研究领域，具体涉及一种基于辅助树的多层快速多极子并行网格细剖方法。

背景技术

[0002] 矩量法是计算电磁学中的一种精确算法，求解的是积分方程。因自动满足辐射边界条件，尤其适用于求解开域问题，如散射和辐射问题。矩量法最终矩阵方程的迭代求解可以采用多层快速多极子技术加速。多层快速多极子技术将目标划分为多层的盒子，以分组分层的方式，通过聚集、转移、发散实现矩阵矢量乘。借助于高效的并行计算技术，可以将矩量法的计算规模扩展到上百亿未知量、上万波长。

[0003] 在面积分方程问题的多层快速多极子技术求解过程中，首先对整个目标的外表面使用三角形单元进行离散，此过程称为网格剖分。为便于计算，对三角形网格剖分信息通常使用以下数组进行存储：

[0004] (1)顶点坐标数组xyz(3，maxnode)，maxnode是剖分产生的总顶点数。每个顶点采用3个实数型元素存储xyz坐标。按照单精度4个字节存储，则存储所需总字节数为12×maxnode。

[0005] (2)三角形单元局部编号与全局编号对应数组ipat(3,maxpatch)，maxpatch是剖分产生的三角形单元总数。每个三角形含3个整形元素存储编号值。按照长整型存储，则存储所需总字节数为24×maxpatch。

[0006] (3)棱边的编号与顶点和三角形对应关系iedge(4,maxedge)，maxedge是剖分产生的总边数。每条边4个整型元素存储。按照长整型存储，则存储所需总字节数为32×maxedge。

[0007] (4)棱边的中点坐标，用于辅助构建多层快速多极子树结构middlenode(3,maxedge),每条边三个实数型元素存储。按照单精度4个字节存储，则存储所需总字节数为12×maxedge。

[0008] 当计算目标电尺寸很大，则离散未知量(三角形棱边数目)规模很大，比如百亿规模时，基于矩量法的网格生成将变得非常困难。这主要有以下原因：

[0009] (1)生成困难。按照以上几何信息存储方式估计，即便是采用二进制存储，最终生成的几何文件大小总量也将接近1TB。考虑中间变量数组，单只网格生成所需的内存就接近甚至超过1TB，无论是硬件资源还是计算时间上考虑，串行程序都难以实现。

[0010] (2)读取困难。网格文件生成后，近TB甚至TB以上量级大文件的读取也将成为并行程序瓶颈，采取单一进程连续读取，网络发送到其他进程的方式效率极低。多极子并行时离散是以多极子树结构为对象进行的，初始剖分网格连续的边编号与按照多极子树结构重排序的边编号顺序不同。在多极子最细层以盒子为单位的近相互作用矩阵填充和聚集、发散矩阵填充时，各进程以多层快速多极子重排后的边编号为偏移量进行几何文件读取，文件指针的偏移数组有序，但不连续。不连续的文件指针偏移将带来磁盘磁头的跳动读取，而在现有的硬盘架构下，此种文件读取方式的读取速度将远慢于一整块等量数据读取。即便是采用并行读取，当进程规模高达数千甚至上万，因硬盘磁头数量限制，读取受限。当大量文件读取工作同时进行，特别是在高性能集群上用户众多，其他用户同时进行文件读取的情况下，几何文件读取极慢。

[0011] (3)高效并行困难。多层快速多极子技术中，为保证计算精度，通常需要保持三角形的平均边长在十分之一波长甚至更小，此种剖分尺度可称为计算离散尺度。此外，还需对目标几何外形精确模拟，此种剖分尺度可称为几何模拟尺度。一般计算中设置几何模拟尺度与计算离散尺度为相同值。然而对于实际目标，几何模拟尺度可远大于计算离散尺度而保持良好的几何外形模拟精度。这为生成初始粗网格，进而通过并行细剖生成满足计算要求的大规模细网格的方法实施成功奠定了基础。受多极子中棱边需重排这一操作的限制，网格细剖实施无法采用一般的基于目标几何形状图划分(metis)的分解模式，否则不合理的初始网格离散，将导致大量的通信。极端情况下，进程上离散的网格需要进行彻底的重分配。

[0012] 综上，大规模网格的生成、读取已经成为并行多层快速多极子程序的瓶颈，迫切需要开发新的技术来实现网格的高效并行生成。

发明内容

[0013] 有鉴于此，本发明提供了一种基于辅助树的多层快速多极子并行网格细剖方法，能够实现大规模网格高效地生成和读取。

[0014] 实现本发明的技术方案如下：

[0015] 一种基于辅助树的多层快速多极子并行网格细剖方法，包括以下步骤：

[0016] 步骤一、将复杂结构目标剖分为三角形网格，根据三角形网格的顶点编号和顶点坐标计算三角形中心坐标，基于三角形中心，按照多层快速多极子盒子划分原则，生成网格的辅助树结构；

[0017] 步骤二、根据预设的并行方案和MPI进程数，对网格进行负载分配，确定归属于各进程上的三角形网格编号；

[0018] 步骤三、对分配后的三角形网格编号、顶点编号和棱边编号进行重排序，并输出对应的存储文件；

[0019] 步骤四、并行读取上述步骤三得到的存储文件，对三角形网格进行预设次数的细剖，并对新生成的棱边、顶点和三角形网格重排序，形成完整的细网格信息。

[0020] 进一步地，步骤三中，对三角形网格编号、顶点编号和棱边编号进行重排序具体为：

[0021] 三角形网格编号：首先对离散到进程上不连续的三角形单元编号按照原编号大小顺序压缩重排为自1开始的连续编号，之后按照进程增序连续编号；

[0022] 顶点编号：先将所有进程内的顶点编号，然后将公共顶点按顺序重新编号；

[0023] 棱边编号：对进程内的棱边按原编号压缩重排后按照进程增序编号，然后对所有公共边压缩重排后编号，公共边编号的起始位置在所有进程内边编号之后。

[0024] 进一步地，步骤四中，细剖的剖分尺度m与次数按如下方式确定：

[0025] m≤m0×2n

[0026] 其中，n是预设的细剖次数，m0是最终需要达到的平均剖分边长，一般设置为十分之一波长。

[0027] 进一步地，步骤四中，对新生成的棱边、顶点和三角形网格重排序，具体为：

[0028] 三角形网格编号和顶点编号采用所有进程的全局编号；

[0029] 顶点编号的原有顶点编号值不变，新产生的顶点按照从属边编号压缩后局部顺序重新编号。

[0030] 有益效果：

[0031] (1)本发明以粗网格三角形中心坐标构建辅助树，以与计算多极子树相同的离散模式对辅助树进行离散，进而确定粗网格单元的离散模式。此种离散方式可保证形成的网格数据文件按进程连续分块，提高文件并行读取效率，更为重要的是，能很好的预估细剖后多极子树的离散模式，从而最大程度的保证细剖后多极子各矩阵填充所需网格信息的本地性，减少进程间数据重分配带来的额外通信。

[0032] (2)本发明粗网格离散到进程后，通过合理设置单元、边、顶点的编号规则，保证后续细剖过程中各进程间完全独立进行，无需通信，因此具有非常高的并行效率，生成网格耗时几乎可忽略，具有高效性。

[0033] (3)网格并行细剖过程中的中间数组在网格细剖操作完成后完全释放，且这些中间数组存储所需内存在何种情况下都远小于后续计算过程的峰值内存，因此不会导致整个并行程序内存需求的增加。

附图说明

[0034] 图1粗网格构建的辅助树离散到进程示意图。

[0035] 图2粗网格几何信息重编号规则示意图。

[0036] 图3三角形单元内顶点与边编号对应关系。

[0037] 图4不同粗网格几何文件并行读取示意图。

[0038] 图5细剖后顶点、棱边、三角形单元编号规则示意图；其中，(a)为细剖后生成三角形单元编号规则示意图，(b)为细剖后生成棱边编号规则示意图，(c)为细剖后生成顶点编号规则示意图。

[0039] 图6细剖后新生顶点坐标计算方法示意图。

[0040] 图7边信息数组iedge组成示意图。

[0041] 图8进程间通信完善iedge数组信息示意图。

[0042] 图9舰船模型几何结构示意图。

[0043] 图10直径2400波长的金属球目标矩阵填充过程中，本地几何信息百分占比对应的进程数统计示意图。

[0044] 图11尺寸6000波长的舰船目标矩阵填充过程中，本地几何信息百分占比对应的进程数统计示意图。

[0045] 图12直径2400波长的金属球VV极化双站RCS与解析解对比示意图。

[0046] 图13舰船目标在10GHz频率下VV、VH极化双站RCS示意图。

[0047] 图14本发明方法流程示意图。

具体实施方式

[0048] 下面结合附图并举实施例，对本发明进行详细描述。

[0049] 本发明提供了一种基于辅助树的多层快速多极子并行网格细剖方法，该技术首先在专业软件生成的初始剖分网格基础上经2-3次细剖形成粗网格，以粗网格三角形单元中心作为未知量空间位置，构建辅助树。对此辅助树采取与计算所用多极子树结构相同的离散模式、进程数进行分层离散，并根据按盒子并行首层的离散模式，经遍历到最细层建立三角形单元与计算所用进程间的分布映射。在此基础上，对三角形单元、边、顶点编号重排，形成分区连续的粗网格文件。其次，并行读取粗网格文件，各进程上对粗网格进行均匀一致性细剖，并对新生成的边、点、三角形按规则编号，形成完整细网格信息。最后，多极子树结构离散后各进程最细层远相互作用聚集/发散矩阵填充、近相互作用矩阵填充时所需的几何信息，通过MPI进程间网络通信完善。数值算例表明，此基于辅助树的多层快速多极子网格并行细剖方法可最大限度的保持数据的本地性，缩减通信数据量，因此具有很高的并行效率。

[0050] 本发明的具体步骤如下：

[0051] 一种基于辅助树的多层快速多极子并行网格细剖方法，如图14所示，包括以下步骤：

[0052] 步骤一、将复杂结构目标剖分为三角形网格，根据三角形网格的顶点编号和顶点坐标计算三角形中心坐标，基于三角形中心，按照多层快速多极子盒子划分原则，生成网格的辅助树结构；

[0053] 步骤1.1：形成初始网格剖分。使用CATIA等商业CAD软件对复杂结构目标进行剖分生成初始网格，将目标剖分成多个面元，并获取每个面元编号、顶点位置坐标。

[0054] 步骤1.2：初始网格细剖生成粗网格。一般商业软件在普通PC机上产生的网格剖分很难达到上千万，为了使辅助树更好的逼近最终细剖后网格建立的多极子树结构，需要通过n1(通常取值2或3)次细剖生成网格量较大的粗网格。

[0055] 步骤1.3：读取粗网格三角形单元顶点编号、顶点坐标，计算三角形中心坐标，三角形中心坐标的计算公式为：

[0056] r＝(r1+r2+r3)/3.0

[0057] 基于三角形中心点，按照多层快速多极子盒子划分思路，生成粗网格的辅助树结构。在此采用自上向下的划分方式，也即构建包围目标的最小尺寸盒子，此为第0层。之后，逐层将盒子在空间三个维度上二分，直到最底层盒子的尺寸小于给定值。实际计算时，通常设置最底层盒子尺寸大小为0.3波长。在辅助树构建时，此值可按比例扩展，也即设置为[0058]

[0059] 其中，n2为由粗网格细剖到最终细网格的细剖次数，设n是总计划细剖次数，于是有：

[0060] n＝n1+n2

[0061] 每个盒子采用空间交织的莫顿码编号，各层盒子按照莫顿码排序顺序存储。本层盒子所对应父层盒子的编号可通过如下计算获得：

[0062] Mi-1＝Mi＞＞3

[0063] 自最底层盒子开始，向上逐层遍历，可确定各层盒子的非空盒子数，进而构建出多层快速多极子八叉树结构。

[0064] 步骤二、根据预设的并行方案和MPI进程数，对网格进行负载分配，确定归属于各进程上的三角形网格编号；

[0065] 如图1所示，根据预先设定的并行方案和计算使用的MPI进程数，进行负载分配，确定辅助树结构的并行层起始层数Lb。之后自该层递归遍历到最细层，确定归属于各进程上的三角形单元编号序列。

[0066] 步骤三、对分配后的三角形网格编号、顶点编号和棱边编号进行重排序，并输出对应的存储文件；

[0067] 几何信息重编号输出至文件。对原网格剖分软件输出的三角形编号、顶点编号、棱边的编号进行重排序。具体的：

[0068] (1)单元编号：每个三角形单元仅能离散到一个进程上，因此首先对离散到进程上不连续的三角形单元编号按照原编号大小顺序压缩重排为自1开始的连续编号，之后按照进程增序连续编号；

[0069] (2)棱边编号：每条边可被两个三角形单元共用，而两个单元可属于不同进程。因此每个进程内的边可分为内部边eI与交界上Γ的共用边eΓ两类。对各进程内边eI的边按原编号压缩重排后按照进程增序编号，最后对所有公共边压缩重排后编号，编号起始位置在所有进程内边编号之后。

[0070] (3)顶点编号：顶点与边编号过程类似，也是按照先所有进程内顶点编号，后公共顶点编号的顺序重新编号。

[0071] 上述过程如图2所示。图中 分别代表i进程上内部边总数，三角形单元总数，内部顶点总数， 和 分别代表公共顶点总数和公共边总数。初始网格的规模通常比较小，因此以上过程可作为前处理过程，串行程序执行即可。最终的输出文件包括：

[0072] Target.xyz文件，记录重排序后顶点对应的坐标值；

[0073] Target.ipat文件，记录每个三角形单元对应的三个顶点全局编号，三条边全局编号。其中三角形内边与顶点编号的对应关系如图3所示；

[0074] Target.shinfo文件，记录公共信息，包括每个进程内局部边、点、三角形数，公共边所属两进程的进程号，规模很小；

[0075] 步骤四、并行读取上述步骤三得到的存储文件，对三角形网格进行预设次数的细剖，并对新生成的棱边、顶点和三角形网格重排序，形成完整的细网格信息。

[0076] 步骤4.1：采取不同方式并行读取步骤五产生的几何文件。其中，[0077] (1)Target.shinfo文件：由0进程读取，广播发送到所有进程；

[0078] (2)Target.ipat文件，各进程独立读取离散到本进程的三角形单元所对应顶点、三条边编号信息，存储至ipat数组；

[0079] (3)Target.xyz文件，各进程独立读取本进程内部顶点所对应行的信息，后0进程读取公共顶点信息，并将广播其发送到所有的进程上，存储为xyz数组。

[0080] 上述文件读取过程如图4所示。从图中我们可以看出，对于规模比较大的Target.xyz和Target.ipat文件，各进程都是一整块的读取从属本地的部分。公共交界顶点部分由0进程读取，依然是一整块数据且公共交界顶点仅为一维线分布，增加速度远小于二维面，因此规模几乎可忽略。

[0081] 步骤4.2：三角形顶点编号局部化。对本进程上的三角形单元所有顶点编号(ipat数组前三列)进行快速排序压缩，其中内部顶点在前，公共顶点在后，压缩后的编号自1开始且连续。根据压缩后的编号与原编号的对应关系，从共有节点坐标数组中提取对应点坐标，并与本进程上内部顶点坐标整合存储。

[0082] 步骤4.3：对网格单元进行预定次数的细剖。剖分所用的剖分尺度m与次数按如下方式确定：

[0083] m≤m0×2n

[0084] 其中n是总的细剖次数，m0是最终需要达到的平均剖分边长，一般设置为十分之一波长。

[0085] 为方便处理公共边信息，采取三角形单元一分为四的均匀一致性细剖。每一次细剖都需分别按照不同规则对对此次细剖产生的三角形单元、边、顶点重新编号，生成顶点坐标信息，存储用于下一次细剖。如图5所示。具体的：

[0086] (1)顶点：进程内局部编号。细剖过程中，每一个原三角形三条边中点上生成新的顶点。顶点的编号规则如图5(a)所示。原有顶点编号值不变，新产生的顶点按照从属边编号压缩后局部顺序 重新编号，起始值为本进程本次细剖前顶点数 加1。细剖后新生顶点坐标的计算方法如图6所示。首先自数组ipat前三列中取第i个三角形对应的三个顶点的进程内局部编号值 之后，根据局部编号自xyz数组中取对应点坐标。最后，由点坐标计算边中点坐标，从而得到对应新生成顶点的坐标。边中点坐标的计算公式为：

[0087] r＝(r1+r2)/2.0

[0088] (2)棱边：采用全局编号，如图5(b)所示。设本次细剖单元的三条边全局编号分别为ei,1,ei,2,ei,3。则此单元每进行一次一分四的均匀细剖，将在单元内部产生三条边，原三角形单元三条边将一分二后形成六条边。单元内三条边的编号按单元全局编号顺序依次编号。边细剖产生的两条边则按照边的全局编号顺序依次编号。需要注意的是，为保证两个单元内对同一条边编号的一致性，定义边上新生成的两条边中，靠近全局编号较小顶点的一条先编号，靠近较大顶点的次之。因边被两个三角形共用，采用此种全局编号方式可保证不同进程上公共边编号的一致性。

[0089] (3)单元：采用所有进程的全局编号。设本次细剖单元i的全局编号为pi，则新产生的四个单元局部编号与全局编号的对应关系如图5(c)所示。

[0090] 步骤4.4：细剖完成后，由ipat数组，生成存储最终棱边信息的iedge数组。每一行共有4个元素，其中前两位为边所对应的顶点，后两位为边所对应的三角形，如图7所示。特别需要指出的是，在此，顶点编号为本进程上的局部编号，而边的编号为所有进程上边的全局编号。

[0091] 步骤4.5：进程间通信，完善iedge数组。步骤九完成后，对于交界面上的公共边，其后两位必然有某个值未填充，需要通过对应进程的通信完善。对于细剖后产生的所有公共边，总是可以通过步骤八所述边的编号规则进行查找。结合Target.sinfo文件中存储的初始公共边对应进程信息，可以很便捷的确定细剖后生成公共边所对应两进程信息。为保证细剖后的边与进程编号的一一对应关系，方便后续矩阵填充时确定边信息收、发进程信息，每条公共边中进程编号大的进程将发送其上此边iedge对应三角形信息至进程编号小的进程上存储，并将本进程上对应信息删除。具体过程如图8所示。

[0092] 步骤4.6：数据重排。因顶点、三角形、棱边在各进程上分片连续存储，可方便的建立点、面、边的进程内局部编号与按进程号顺序全局编号间的映射关系。从而为后续通过全局编号确立顶点、边、面诸元素收发进程奠定基础。

[0093] 经以上过程，完成并行网格细剖过程，生成完整细剖网格信息。

[0094] 之后，按照一般并行多层快速多极子程序执行流程进行计算。在并行多层快速多极子近相互作用矩阵填充、最细层聚集发散矩阵填充等需要使用几何信息过程时，非本地所需的几何信息可通过确定进程、局部编号的方式，经网络通信获得。因几何信息的离散时使用了辅助树结构，因此可以保证绝大部分的几何信息都处于本进程。仅有少部分信息需要通信完成，从而保证程序效率。

[0095] 实施例：

[0096] 本实施例以直径为2400波长的金属球和如图9所示长约6000波长的舰船目标模型为目标，采用本发明提出的并行网格细剖方法，调用960个MPI进程进行细剖后调用三元并行多层快速多极子程序计算，以测试生成网格的正确性。计算平台为中科院网络中心“元”超级计算平台，每个节点256G内存，2个Intel-E5-2680V3处理器，24CPU核。

[0097] 对两个目标细剖生成网格信息如表1所示。其中初始网格是由CATIA剖分产生，粗网格为构建辅助树所用的由初始网格经3次细剖产生，最终网格为满足计算所需尺度的网格剖分。从中可以看出，并行细剖到几十亿网格所需时间仅为10余秒，相比于数小时的整体计算时间基本可以忽略。

[0098] 表1金属球和舰船目标并行细剖信息统计表

[0099]

[0100] 统计上述两个不同目标，在远相互作用最细层盒子聚集、发散矩阵填充过程中，本地存储几何信息数据量占总共所需数据量的百分比情况。以边信息为例，图10是金属球目标时，对所有960个进程进行统计，拥有不同本地数据百分占比进程数的情况。从图中可以看出，绝大多数进程上计算所需几何信息文件的本地存储数据占比都在95％以上，表明这些进程在填充聚集发散矩阵时，仅有5％的数据需要通过通信获得，本地数据约为需通信发送数据的20倍，大一个量级以上。同理，对于舰船目标，如图11所示，对所有960个进程进行统计，拥有不同本地数据百分占比进程数的情况，可以得到类似结论。因此，辅助树的构建，可以很好的近似多极子八叉树结构，保证离散存储的良好预分配，达到缩减额外通信，提高效率的目的。

[0101] 金属球目标多层快速多极子程序计算输出的双站RCS结果(PMLFMA)与解析解(Mie)对比如图12所示，雷达入射波的入射角度为θ＝0°， 观察平面为xoy平面，其中θ＝0°为后向散射。经统计，两者之间的均方根误差(RMS)为0.7dB,说明了提出的基于辅助树的多层快速多极子网格并行细剖方法的正确性、有效性。对舰船目标的计算结果如图13所示，入射角度为θ＝70°， 观察角度为θ＝70°， 说明了此并行网格细剖方法对复杂目标同样有效，具有通用性。

[0102] 本发明实施例所展示的直径2400波长金属球、6000波长舰船模型经两次细剖生成计算所需数十亿网格，计算结果正确有效，证明了本发明并行方法的实用性，可实现大规模网格并行细剖生成。

[0103] 综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。