摘　要　提出了一种基于网格生成递归法的并行区域划分算法，该算法依据网格生成代价的估算分析，采用迭代分解法对区域进行并行划分.在曙光1000A 系统上的运行结果表明，该网格算法的效率和加速比均优于串行递归算法.
    关键词　有限元网格；并行区域划分算法；网格生成代价；迭代分解法

　　基于网络生成递归法^［1～3］，本文提出了一种并行区域划分算法，该算法满足以下四个基本原则：a. 任务平衡性原则.能生成平衡的子区域集，即在各子区域中生成网格的时间大致相等.b. 边界最简原则.子区域的边界结构简单，边界处理所需时间短，处理器间消息传递的费用低.c. 网格均匀原则.并行生成的最终网格形状均匀，无奇异单元.d. 区域划分代价最小原则.区域划分算法本身的代价尽可能小.

1、基本思想及相关算法
　　在网格生成递归法中，如果每个子区域都包含相同的单元数，就比较容易实现任务平衡.因此，首先按照单元数估算待处理区域的网格生成代价，然后根据当前参与并行处理的处理器数N对区域进行分解，并对分解所得子区域进行边界处理，最终获得相互之间既平衡又独立的N个并行子任务.
1.1　网格生成代价的估算算法
　　网格生成代价与分布于待处理区域中的单元数目紧密相关，而单元数目是由该区域的总面积S和区域内单元分布密度决定的.估算公式如下：

G=S/S_tri,　（1）

S_tri=［L²/(2M²)］sin60°,　　（2）

M=Σ2l_i/(d_i1+d_i2) (i=1, 2, …， n),　　　（3）

式中，G表示三维平面区域 (含有n条边) 的网格生成代价；S表示该区域的面积；Stri 是对区域内单个三角形单元平均面积的粗略估算值；L表示区域的周长；M表示区域边界上总的离散段数；l_i表示第i条边的长度;d_i1和d_i2 分别表示第i条边前端点和后端点处的单元边长.
　　据以上公式，网格生成代价的估算算法步骤如下.
　　步骤 1　依次计算区域各边界的长度和离散段数,根据式(2)和(3)计算S_tri；
　　步骤 2　应用三角累加算法计算S,根据式(1)求得G.三角累加算法步骤如下：a. 置变量S为0；b. 选取区域的任意一个顶点作为V₁,按顺时针 (或逆时针) 方向取与V₁ 依次相邻的两个点为V₂和V₃;c. 计算由V₁, V₂ 和 V₃ 组成的三角形的面积，累加入变量S;d. 以V₃为新的V₂，按顺时针 (或逆时针) 方向取与V₂相邻的第一个点为新的V₃，若V₃=V₁，则算法停止；否则转 c.
1.2　区域划分算法
　　估算出整个区域的网格生成代价后，区域划分算法的任务就是：寻找N-1条分割线，将区域分割为N个子区域，使得各个子区域的网格生成代价大致相等.由于无法精确确定子区域中的单元数目，因此允许各子区域的网格生成代价在G/N附近有一个δ误差，亦即所产生的子区域的网格生成代价都属于［G/N-δ, G/N+δ］区间.
　　根据原则 b，任何一条分割线的两端点都定位在区域的边界上，而不落在区域中.在图 1 （a） 中，任何一个子域都只通过一条完整的分割线与其他子域相邻；而在图1（b） 中，子域3与子域1和2之间的相邻关系不便于边界处理.

　　　　　　(a)　　　　　　　　　(b)
图 1　区域划分效果图

　　根据原则c，多条分割线不交于同一点，因为这种分割在N较大的情况下会导致极点现象——在多条分割线的交点处出现小内角的长薄单元，影响网格形状.
　　本文采用迭代分解法划分区域，步骤如下.
　　步骤 1　记待分解的区域为 D;
　　步骤 2　循环执行以下步骤N-2次:a. 在D中找一条分割线将区域分为两个子域D₁和D₂，使D₁的网格生成代价约等于G/N;b. 将D₁从D中去掉，以D₂为新的D，转步骤 a;
　　步骤 3　在D中找一条分割线将区域分为两个子域D₁和D₂，使D₁和D₂的网格生成代价均属于 ［G/N-δ,G/N+δ］.
　　本算法中，寻找分割线的方法是：先在区域中找到一个合适的顶点作为分割线的前端点X，然后从这个顶点出发，沿顺时针（或逆时针)方向按对半查找法在区域的边界上寻找分割线的后端点Y.对半查找算法步骤如下：
　　步骤 1　从X点出发，按顺时针（或逆时针）方向取与X依次相邻的两个点为R和Y；
　　步骤 2　计算分割线XY左侧 （或右侧） 子域的网格生成代价g;
　　步骤 3　若g∈［G/N-δ,G/N+δ］，则算法停止；若g>G/N+δ，则令T=Y，以T与R的中点为新的Y，转步骤2；若g<G/N-δ,则令R=Y.若Y是原区域的顶点，则沿顺时针 （或逆时针） 方向取与Y相邻的第一个点为新的Y，否则，以R与T的中点为新的Y，转步骤2.
　　根据以下原则确定分割线的前端点X：对第一次分割，选取区域中最大顶角所在顶点为X；对第 i （i=2, 3, …, N-1） 次分割，选取分割线两端点中较大顶角所在的顶点为X.根据这一策略，每次分割都选取顶角较大的点作为前端点，一定程度地避免了极点的形成，但还不能完全避免极点的产生.对此，引入分割线共点控制机制：首先定义参数R，由它指明能交于一点的分割线的最大数目,然后，每次分割都用变量r记录X处的共点率，一旦发现r>R，就变换X.
　　图 2 是对区域ABCDEF作6次划分的示意图.在图2（a）中，R取值为2，寻找第5条分割线时，因为∠CHI比∠HIE大,X原应定位在H点，但由于已有两条分割线经过H点，因此如果这次还选取H作为前端点， 那么共点率r增为3,超过预定义的R值,于是，改选I点为X；而在图2(b)中，因为R取值为3，所以，H仍可作为第5条分割线的前端点.

　　　　　　(a) R=2　　　　　　　　　(b) R=3
图 2　共点控制

　　该区域划分算法的特点是：并行任务的负载平衡程度由δ调节；所有分割点都定位在区域的边界上，因为每次分割剩下的子域只有一条边是在分割过程中新加入的分割线，无论选取这条边的哪个端点作为X，X和Y总是落在区域的边界上，有效地满足了边界最简原则,；通过R灵活控制分割线共点率，消除极点，保证网格均匀度.

2、并行实现
　　本文研究的并行区域划分算法已在曙光1000A系统上实现.根据曙光1000A系统目前的条件和特点，选取PVM作为并行编程环境，分三个步骤实现有限元网格的并行生成：首先，由node₁上的master任务程序检测当前PVM中的结点数N，运用区域划分算法分解待处理的区域为N个子区域，将各子区域分配到各node上；然后，并行执行各处理器上的slave任务程序，在各子区域中生成均匀网格；最后，由node₁上的master任务程序收集、组装在各处理器上生成的子区域网格.这里，采用了动态负载平衡策略——farm模式，但有一点与通常情况下不同：为了提高处理器利用率，node₁在子区域分配过程中也分得一个子区域，所以，在数据分配之后，子区域网格传回之前，主处理器也不空闲.
　　由于在并行划分阶段已经对边界进行了预处理，因此，一旦子区域分配给各个处理器，所有处理器将同时独立地在自己的数据子块上根据预定义的单元分布密度生成网格.经过大致相等的时间段后，它们将完成各自的工作，并将生成的子区域网格传送给master任务程序.整个并行算法的通信很小，仅包括从node₁向其他各node广播子区域数据，以及各node向node₁回送生成的子区域网格单元.

3、实验结果
　　表 1 以实例说明上述并行算法在曙光1000A系统上的运行效果 （表中括号内外的数据分别对应例 1 和例 2）.在例 1 中，δ＝G/(5N(N-1))，G=15　987 个,实际生成的单元总数为12　531个；例 2 中，δ＝G/(5N(N-1)),G＝33　452 个，实际生成的单元总数为28　769个.由于区域划分和组装数据所耗费的时间是μs级的，因此将并行算法的执行时间 (t) 分为两部分：子区域网格生成时间 (t₁) 和通信时间 (t₂).t₁ 指各node在其分配到的子区域中生成网格的时间；t₂包括node₁向其他各node发送子区域数据的时间，以及各node向node₁回送子区域网格数据的时间.表 1 给出了使用2,3,4,5,6个node进行并行网格生成的时间开销及加速比 (ν).使用一个处理器时，既不需要做区域划分，也不用通信，因此，在区域中生成网格的时间即为全部运行时间.

表 1　并行算法执行时间及加速比

　　表 2 给出了例1中并行算法在六个node上运行时各node上生成的网格单元总数的估算值 (n) 和实际值 (n?)，以及在子区域中生成网格的时间 (t₃).由表 2 可见，各子区域单元总数的估算值与实际值的比例分布一致，各node之间网格生成时间差别小于5　%，达到了很好的负载平衡.值得注意的是，网格生成时间并不总是与网格单元数目成正比，其原因在于网格生成的时间开销不仅受子区域内网格单元数目的影响，而且与网格单元在子区域中的分布情况有关.

表 2　负载分配表

参考文献
　1　Facello M A. Implementation of a Randomized Algorithm for Delaunery and Regular Triangulations in Three Dimensions. Computer Aided Geometric Design, 1995 (12): 349～370
　2　Chan C T, Anastasion K. An Automatic Tetrahedral Mesh Generation Scheme by the Advancing Front Method. Communications in Numerical Methods in Engineering, 1997, 　13: 33～46
　3　Zeng Yong, Cheng Gengdong. Knowledge-Based Free Mesh Generation of Quadrilateral Elements in Two-Dimensional Domains. Microcomputers in Civil Engineering, 1993 (8): 259～270