随着先进制造与智能制造的高速发展，几何设计与仿真的无缝数据融合已成为CAD/CAE领域发展中亟需解决的关键瓶颈问题。在CAE领域的产品仿真分析阶段，无论是经典的有限元分析方法，还是新兴的等几何分析方法，计算域的网格剖分质量均是影响其仿真结果精度和计算效率的重要因素。目前，网格生成已经成为计算机辅助工程一个重要的研究领域^[1]，在工程与图形学中具有非常关键的应用价值。

三角网格与四面体网格是有限元分析中最为常用的计算域几何表示形式，也是网格生成主要的研究对象。同三角网格相比，四边网格的连接关系更加规则，在整体结构和单元性质方面具有天然的优势^[2]，而且基于四边形网格的有限元求解能达到更高的精度；另一方面，等几何分析方法的提出，将计算域的参数化问题推向研究前沿，而高质量的四边形网格剖分，是计算域参数化问题的关键几何基础^[3]。高质量四边形网格的生成，已引起人们越来越多的关注。但由于其几何形状和拓扑连接关系的强约束性，如何自动、鲁棒、高质量地生成四边形网格，是困扰网格生成领域的难题之一。四边网格生成方法的研究主要有两大学派：其中一部分是从有限元仿真模拟的角度研究平面区域的四边网格化问题；另外一部分则是从动漫几何造型的视角研究3维网格曲面上的四边网格化构造问题，主要包括参数化映射方法^[4]，Morse-Smale复形方法^[5]，方向场方法等。本文所研究的问题属于第一部分的研究工作，即研究平面区域的四边形网格生成，目前在商业软件中所采用的主流方法主要包括三角网格转换法^[6]、铺砌法^[7]、模板法^[8-10]、中轴法^[11]等。三角网格转换方法通常会受限于三角网格的顶点分布和顶点之间的连接关系，难以生成奇异点较少的高质量四边形网格^[6]；基于中轴分解的四边网格生成方法对于边界变化较为敏感，算法整体上不够鲁棒且难以自动化^[11]; 针对复杂平面区域，铺砌法^[7]可以产生高质量四边网格，但由于在不同方向上推进的单元相遇时可能会发生自交现象，因而无法保证所生成网格的有效性，另一方面，对于复杂边界，该方法会在内部生成较多的奇异点。

针对目前平面区域四边形网格化方法的上述局限性，本文面向CAD/CAE数据无缝集成的迫切需求，提出了一种基于边界简化与多目标优化的高质量四边形网格生成新框架，与现有方法相比，本文方法具有较好的鲁棒性，能自动生成有效的全四边形网格，另外，本方法所生成的四边形网格不仅具有较少的网格顶点和单元总数目, 而且奇异点数目可减少70%80%，具有良好的光滑性、均匀性和正交性，非常适用于工程分析和动画仿真。所给出的若干比较实例验证了本文方法的有效性。

本文所研究的问题可描述如下：给定2维平面中的离散折线边界，如何构造以该折线为边界的四边网格，使其网格质量满足工程需求。

针对高质量四边形网格生成这一问题，拟采取的算法可概括为图 1中的5步。首先针对亏格非零的平面区域，提出了一种将多连通区域转化为单连通区域的方法，可生成高质量的插入边界；然后提出了“可简化角度”和“可简化面积比率”两个阈值概念，从顶点夹角和顶点三角形面积入手，将给定的多边形边界简化为粗糙多边形；然后，对边界简化得到的粗糙多边形进行子域分解，并根据奇异点较少的原则，确定每个子域内的网格顶点连接信息；最后通过多目标非线性优化技术确定网格内部顶点的几何位置，可生成具有良好光滑性、均匀性和正交性的高质量四边形网格。下面将分别对每一步骤的细节详细阐述。

一条封闭边界可表示为$\mathit{\boldsymbol{M}} = \left( {\mathit{\boldsymbol{V}}, \mathit{\boldsymbol{E}}} \right)$，其中$\mathit{\boldsymbol{V}} = \left( {{v_0}, {v_1}, \cdots, {v_{n - 1}}} \right)$表示所有的顶点，$\mathit{\boldsymbol{E}}$表示边集，其集合的元素为$\left( {{v_i}, {v_j}} \right)$，其中${v_i}, {v_j} \in \mathit{\boldsymbol{V}}$，$j = \left( {i + 1} \right){\rm{ mod}}\;\mathit{n}$。对于一个多连通域，可表示为一条外边界${\mathit{\boldsymbol{M}}_1}$和$n-1$条内边界${\mathit{\boldsymbol{M}}_2}, \cdots, {\mathit{\boldsymbol{M}}_n}$。对于此类边界，需要用一个二元组$\mathit{\boldsymbol{M}}$来等价代替${\mathit{\boldsymbol{M}}_1}, {\mathit{\boldsymbol{M}}_2}, \cdots, {\mathit{\boldsymbol{M}}_n}$。

以只有一条内边界的多连通域为例，即对于输入的多连通区域总是由两条边界${\mathit{\boldsymbol{M}}_1} = \left( {{\mathit{\boldsymbol{V}}_1}, {\mathit{\boldsymbol{E}}_1}} \right)$，${\mathit{\boldsymbol{M}}_2} = \left( {{\mathit{\boldsymbol{V}}_2}, {\mathit{\boldsymbol{E}}_2}} \right)$组成(如图 2(a)所示)。首先在顶点集${\mathit{\boldsymbol{V}}_1}$和${\mathit{\boldsymbol{V}}_2}$中找到顶点${v_i}, {u_j}\left( {{v_i} \in {\mathit{\boldsymbol{V}}_1}, {u_j} \in {\mathit{\boldsymbol{V}}_2}} \right)$，使其满足该两点间的距离$l$是所有内边界顶点到外边界顶点的距离最小，即

$ l\left( {{v_i}, {u_j}} \right) = {\rm{min}}\left( {l\left( {{v_i}, {u_j}} \right)} \right), {v_i} \in {\mathit{\boldsymbol{V}}_1}, {u_j} \in {\mathit{\boldsymbol{V}}_2} $

在得到${v_i}, {u_j}$之后，即可由顶点集${\mathit{\boldsymbol{V}}_1}, {\mathit{\boldsymbol{V}}_2}$

$ \begin{array}{l} {\mathit{\boldsymbol{V}}_1} = \left\{ {{v_0}, {v_1}, \cdots, {v_{n1 - 1}}} \right\}\\ {\mathit{\boldsymbol{V}}_2} = \left\{ {{u_0}, {u_1}, \cdots, {u_{n2 - 1}}} \right\} \end{array} $

得出顶点集$\mathit{\boldsymbol{V}}$的序列

$ \mathit{\boldsymbol{V}} = \left\{ {{v_0}, {v_1}, \cdots, {v_i}, {u_j}, \cdots, {u_j}, {v_i}, \cdots, {v_{n - 1}}} \right\} $

式中，${u_j}, \cdots, {u_j}$为顺时针方向。

可以看出，边$\left( {{v_i}, {u_j}} \right)$和$\left( {{u_j}, {v_i}} \right)$是一对新加入的位置重叠的边，集合$\mathit{\boldsymbol{V}}$的顶点个数则为$n = {n_1} + {n_2} + 2$。为使新的单连通边界$\mathit{\boldsymbol{M}}$的顶点集$\mathit{\boldsymbol{V}}$中顶点间的疏密程度保持均匀，则还需要对边$\left( {{v_i}, {u_j}} \right)$和$\left( {{u_j}, {v_i}} \right)$进行均匀插值，即在顶点${{v_i}}$到${{u_j}}$和${{u_j}}$到${{v_i}}$中加入新的顶点${{z_k}}$。加入的顶点个数$m$则为${{v_i}}$到${{u_j}}$的距离与边集合${\mathit{\boldsymbol{E}}_1}$和${\mathit{\boldsymbol{E}}_2}$中所有边的平均长度之比值，即

$ \left| {\frac{\begin{array}{l} \;\;\;\;\;\;\;\;\;\;m \approx \\ \left( {{n_1} + {n_2}} \right) \times l\left( {{v_i}, {u_j}} \right) \end{array}}{{\left[{\sum\limits_{i = 0, j = (i + 1){\rm{mod}}{\mathit{n}_1}}^{i = {n_1}-1} {l\left( {{v_i}, {v_j}} \right)} + \sum\limits_{i = 0, j = (i + 1){\rm{mod}}{\mathit{n}_2}}^{i = {n_2}-1} {l\left( {{u_i}, {u_j}} \right)} } \right]}}} \right| $

然后再将这些顶点均匀插值到顶点${v_i}$到顶点${u_j}$和顶点${u_j}$到顶点${v_i}$之间。则变换成单连通的边界$\mathit{\boldsymbol{M}}$的顶点集$\mathit{\boldsymbol{v}}$的序列为

$ \mathit{\boldsymbol{V}} = \left. {\left\{ {\begin{array}{*{20}{l}} {{v_0},{v_1}, \cdots ,{v_i},{z_0}, \cdots ,{z_{m - 1}},{u_j}, \cdots ,}\\ {{u_j},{z_{m - 1}}, \cdots ,{z_0},{v_i}, \cdots ,{v_{n - 1}}} \end{array}} \right.} \right\}\ $

若给定的边界含有较多的顶点，则会在子域分割的可视点链计算过程中带来大量的计算，而其中多数相邻两顶点的计算结果是不必要的。本文提出在子域划分前首先进行边界简化这一思路，以合理地减少顶点数量，降低计算量；其次，对于简化所得出的几何特征点也为确定子域边界的角点提供了额外信息，从而不必计算边界上每一个顶点。

假设多边形边界用一个二元组$\mathit{\boldsymbol{M}} = \left( {\mathit{\boldsymbol{V}}, \mathit{\boldsymbol{E}}} \right)$表示，其中$\mathit{\boldsymbol{V}} = \left( {{v_0}, {v_1}, \cdots, {v_{n - 1}}} \right)$表示所有的顶点，$\mathit{\boldsymbol{E}}$表示边集，其集合的元素为二元组数对(${v_i}$, ${v_j}$)，其中${v_i}, {v_j} \in \mathit{\boldsymbol{V, }}j = \left( {i + 1} \right){\rm{mod}}\;\mathit{n}$。简化的目标是在尽量保证多边形边界几何形状的前提下，尽可能多的减少边界顶点的个数，最后得到简化后的边界${\mathit{\boldsymbol{M}}_s} = \left( {{\mathit{\boldsymbol{V}}_s}, {\mathit{\boldsymbol{E}}_s}} \right)$，其中${\mathit{\boldsymbol{V}}_s} \subseteq \mathit{\boldsymbol{V}}$。

3维空间网格模型的简化技术已非常成熟。针对2维多边形的边界化简问题，从顶点夹角和顶点三角形面积入手，本文对文献[12]中的简化方法进行了改进，提出了“可简化角度”和“可简化面积比率”两个阈值，并提出两个简化标准用于控制边界化简的程度。相关概念定义如下：

定义1  顶点夹角$\alpha $：在多边形边界上顶点${v_i}$和其相邻的两个顶点${v_{i - 1}}$、${v_{i + 1}}$构成的以${v_i}$为顶点的夹角$\alpha$的大小，其中$0 \le a \le \pi $，如图 3(a)所示。

定义2  顶点三角形面积$S$：在多边形边界上某一个顶点${v_i}$和其相邻的两个顶点${v_{i - 1}}$、${v_{i + 1}}$构成的三角形面积$S$，如图 3(b)所示。

在本文中，基于顶点夹角和顶点三角形面积，通过“可简化角度”和“可简化面积比率”两个阈值来控制简化程度：1)当顶点夹角小于可简化角度$A\left( {0 \le A \le \pi } \right)$时，则该点不能被简化；2)当顶点三角形面积和整个多边形的面积的比率大于可简化面积比率$\gamma \left( {0 \le \gamma \le 1.0} \right)$时，则该点不能被简化。本文所采用的边界化简算法可总结如下：

1) 令${\mathit{\boldsymbol{M}}_s} = \mathit{\boldsymbol{M}}$，即初始简化边界与原边界相同。

2) 对于每一个顶点${v_i}$∈${\mathit{\boldsymbol{V}}_s}$，计算其相应的顶点夹角${\alpha _i}$和顶点三角形面积${S_i}$。

3) 将简化边界中所有顶点的顶点夹角与可简化角度比较，以及将该顶点三角形面积与简化边界所构成的多边形面积的比率与可简化面积比率进行比较。找出所有可简化顶点的中顶点夹角最大的顶点${v'}$。如果没有找到，则简化过程结束，否则执行第4)步。

4) 将${v'}$从${\mathit{\boldsymbol{V}}_s}$集合中去除，继续执行步骤2)。

由边界化简算法可以看出，边界简化程度可通过两个阈值即可简化角度$A$和可简化面积比率$\gamma $来进行控制。本文提出以下两个条件来控制简化程度:

1) 对于简化结束后的边界${\mathit{\boldsymbol{M}}_s} $，每一条边$e_i$∈${\mathit{\boldsymbol{E}}_s}$，该边所连接的两个顶点${v_i}$, ${v_j}$在其原边界$\mathit{\boldsymbol{M}}$中，该两点中间可能会有一系列被简化掉的顶点${z_k}, 0 \le k \le \left( {j - i} \right){\rm{mod}}\;\mathit{n}$。对于这些顶点其$x$坐标或$y$坐标至少必须有一组满足非递减或非递增条件，即$x$坐标满足：${x_{z0}} \le {x_{z1}} \le \cdots \le {x_{zk}}$或${x_{z0}} \ge {x_{z1}} \ge \cdots \ge {x_{zk}}$或$y$坐标满足${y_{z0}} \le {y_{z1}} \le \cdots \le {y_{zk}}$或${y_{z0}} \ge {y_{z1}} \ge \cdots \ge {y_{zk}}$。

2) 对于简化结束后的边界${\mathit{\boldsymbol{M}}_s} $，每一条边$e_i$∈${\mathit{\boldsymbol{E}}_s}$，该边所连接的两个顶点${v_i}$, ${v_j}$在其原边界$\mathit{\boldsymbol{M}}$中，该两点中间可能会有一系列被简化掉的顶点${z_k}, 0 \le k \le \left( {j - i} \right){\rm{mod}}\;\mathit{n}$。这些点的位置必须在一定范围内。以${v_i}$, ${v_j}$为两定点，一动点$p$到该两定点所连线形成的夹角为一定值$\beta $，则该动点$p$的移动轨迹所形成的区域，即为被简化掉的顶点的分布范围，如图 4所示。

需要说明的是，${v_i}$, ${v_j}$之间被简化掉的顶点${z_k}$在满足条件1)的情况下也必须要满足条件2)。另外，可简化角度$A$和可简化面积比率$\gamma $越接近同时满足条件1)和条件2)的临界值越好。图 2中所给例子的边界简化结果如图 5所示。

本节研究如何将一个多边形区域划分为多个几何形状相对简单的多边形区域，降低后续生成拓扑结构、网格映射等步骤的复杂度，并且能大大减少网格奇异顶点数目，提高拓扑网格的质量。

对于给定多边形边界$\mathit{\boldsymbol{M}} = \left( {\mathit{\boldsymbol{V}}, \mathit{\boldsymbol{E}}} \right)$，将其分解为多个多边形边界${\mathit{\boldsymbol{M}}_1} = \left( {{\mathit{\boldsymbol{V}}_1}, {\mathit{\boldsymbol{E}}_1}} \right), \cdots, {\mathit{\boldsymbol{M}}_n} = \left( {{\mathit{\boldsymbol{V}}_n}, {\mathit{\boldsymbol{E}}_n}} \right)$，满足原多边形的几何形状和分解后的多个多边形的整体几何形状相同即可。

分解原边界的输入即为在3.3节多边形边界化简的输出结果和原边界，即图 2(b)和图 5所示。在进行分解前还需要对多边形边界简化后的边界加以处理，即在子域分解过程中对多边形边界数据进行前处理，处理方法如下^[13]：

对于简化结束后的边界${\mathit{\boldsymbol{M}}_s} $，每一条边$e_i$∈${\mathit{\boldsymbol{E}}_s}$，该边所连接的两个顶点${v_i}$, ${v_j}$在其原边界$\mathit{\boldsymbol{M}}$中，该两点中间可能会有一系列被简化掉的顶点${z_k}, {\rm{ }}0 \le k \le \left( {j - i} \right){\rm{mod}}\;\mathit{n}$。${z_k}$满足多边形边界简化的限制条件1)和限制条件4)。${z_k}$中的点在简化边界${\mathit{\boldsymbol{M}}_s} $已经被完全删除，但此时需要重新把这些点加入到简化边界${\mathit{\boldsymbol{M}}_s} $中以构成新的边界，即子域分解输入边界${\mathit{\boldsymbol{M}}_d} $。本文将这些点分为虚顶点和实顶点两类，定义分别如下：

定义3  虚顶点：在子域分解输入边界中${\mathit{\boldsymbol{M}}_d} $中那些位置坐标被重新计算过后的顶点，也就是在化简原边界过程中删除过的那些顶点。

定义4  实顶点：在子域分解输入边界中${\mathit{\boldsymbol{M}}_d} $中那些拥有真实位置信息的顶点，也就是在化简原边界过程中保留下来的那些顶点。

重新计算虚顶点的位置信息的方法如下，对于简化结束后的边界${\mathit{\boldsymbol{M}}_s} $，取其中一条边$e_i$∈${\mathit{\boldsymbol{E}}_s}$，该边所连接的两个顶点${v_i}$, ${v_j}$在其原边界$\mathit{\boldsymbol{M}}$中，该两点中间可能会有一系列被简化掉的顶点${z_k}$, 连接${v_i}$, ${v_j}$构成线段$s$，分别通过每一个顶点${z_k}$对线段$s$做垂线，相交与线段$s$于${{z'}_k}$点，即用${{z'}_k}$点的位置信息来替换原顶点信息。对于化简边界${\mathit{\boldsymbol{M}}_s} $中的每一条都进行这样的处理，最终得出子域分解输入边界${\mathit{\boldsymbol{M}}_d} $。对于如图 5所示简化边界，则其对应的子域分解输入边界${\mathit{\boldsymbol{M}}_d} $则如图 6(a)所示。

得到子域分解输入边界${\mathit{\boldsymbol{M}}_d} $之后，在分解过程中，就不必考虑多边形边上的每一个顶点，即分解每一条子域的起始点只会从边界中的实顶点中选出，结束的顶点可以为实顶点，也可以为虚顶点。对于是分解出的子域边界，如果结束顶点是虚顶点，则将其标记为实顶点。

对于一条子域分解输入边界${\mathit{\boldsymbol{M}}_d} $，首要问题是如何选择一条子域边界的起始顶点，即如何从边界${\mathit{\boldsymbol{M}}_d} $挑选出一个分解起始点$v$。选取依据主要来源于${\mathit{\boldsymbol{M}}_d} $中所有实顶点中为凹顶点的可视度和其可视顶点链信息。可视度可分为正可视度和逆可视度两个方向，正可视度又可分为正实可视度和正虚可视度，同样逆可视度也将分为逆实可视度和逆虚可视度。具体定义如下：

定义5  正(逆)虚可视顶点链：设顶点$v$为分解输入边界${\mathit{\boldsymbol{M}}_d} $的一个实顶点，从顶点$v$逆(顺)时针方向的相邻点开始，到以顶点$v$逆(顺)时针方向下一个实顶点的下一个顶点起与视点$v$可视的连续顶点序列的最后一个顶点为止，该顶点序列即为正(逆)虚可视顶点链，包含实顶点和虚顶点。

定义6  (正(逆)虚可视度)：即正(逆)虚可视顶点链的顶点个数。

定义7  正(逆)实可视顶点链：即正(逆)虚可视顶点链中包含的实顶点序列。

定义8  正(逆)实可视度：即正(逆)实可视顶点链的顶点个数。

对所有凹顶点中实顶点的正逆虚实可视度和可视顶点链的相关数据计算完毕之后，所选取的起始顶点$v$点的实可视度相同，则比较两实顶点的实可视度相关的虚可视度，选择虚可视度相对较大的一个顶点，若虚可视度再次相同，则两个实顶点是等价的，此时选择顶点序列较小的顶点。

为使整个子域边界所有边的长度尽量保持均匀，则需要对由子域边界中的起始顶点和结束顶点${v'}$所构成的边$l$进行插值，插值的顶点个数$n$为边$l$的长度除以子域边界其他边的平均长度。由于需要保证子域边界的顶点个数为偶数，故当插入顶点个数加上子域边界上顶点个数不是偶数时，则插入顶点个数$n$则需要减去1，若插入顶点个数为0，则$n$需要加1。确定好插入顶点个数后，则需要分别对原边界${\mathit{\boldsymbol{M}}_d} $和分解的子域边界${\mathit{\boldsymbol{M}}'}$中的顶点$v$和${v'}$之间进行坐标位置信息的均匀插值。

如此循环，直到最后一次计算得到的子域分解的边界与子域分解输入边界完全相同，则表示分解结束。该次分解结束后下一次子域分解输入边界${{\mathit{\boldsymbol{M'}}}_\mathit{d}}$在插值前将只会包含两个顶点，即上一次分解的子域边界的起始顶点和结束顶点。对于图 6(a)的子域分解结果如图 6(b)所示。

综上所述，子域分解步骤中的算法流程可总结描述如下(如图 7所示)：

1) 根据原边界$\mathit{\boldsymbol{M}}$和简化边界${\mathit{\boldsymbol{M}}_s} $进行子域分解步骤前处理，得到子域分解输入边界${\mathit{\boldsymbol{M}}_d} $。

2) 计算${\mathit{\boldsymbol{M}}_d} $中所有凹顶点中实顶点的可视度信息。

3) 根据得到的可视度信息，在边界${\mathit{\boldsymbol{M}}_d} $中所有实顶点中选择出本次子域分解的起始顶点$v$。

4) 计算以顶点$v$为起点的子域边界链顶点序列。

5) 确定子域边界顶点链的顶点序列之后，就可得到出子域边界${\mathit{\boldsymbol{M'}}}$。将子域边界中除起始顶点和结束顶点之外的其他顶点从边界${\mathit{\boldsymbol{M}}_d} $中删除，得到下一次分解输入边界${{\mathit{\boldsymbol{M'}}}_\mathit{d}}$。

6) 计算子域边界起点顶点$v$和结束顶点${v'}$所构成的边中需要插入的个数

$ n = \frac{{l\left( {{{e'}_{n - 1}}} \right)}}{{\sum\limits_i^{i = n - 2} {l\left( {{{e'}_i}} \right)} }}, e' \in \mathit{\boldsymbol{E}}' $

式中，$l$为子域边界的长度。

如果$n + n'$为不为偶数，且$n>0$则令$n=n-1$，否则$n=1$。

7) 对子域边界 ${\mathit{\boldsymbol{M'}}}$中的顶点$v$和${\mathit{v'}}$间进行均匀插值，插入的顶点的个数为$n$。

8) 如果此时边界${{\mathit{\boldsymbol{M'}}}_\mathit{d}}$中的顶点个数为2，则子域分解结束。否则对边界${{\mathit{\boldsymbol{M'}}}_\mathit{d}}$中的顶点$v$和${\mathit{v'}}$进行均匀插值，插入的顶点的个数同样为$n$。

9) 令${\mathit{\boldsymbol{M}}_d} $=${{\mathit{\boldsymbol{M'}}}_\mathit{d}}$，继续执行第2)步。

图 7中给出了一个子域分解实例及相应的步骤详解。图 7(a)为简化后的输入边界，其中蓝色顶点为实顶点，黄色顶点为虚顶点。首先计算边界$G$中所有凹实顶点的可视度信息, 并将此信息标记为图中的格式4(22)-$v_3$-4(23), 其含义为：逆实可视度(逆虚可视度)—顶点—正实可视度(正虚可视度)。由此确定$v_3$分解为起始点, 分解方向为逆时针; 图 7(b)为确定第1次分解后得到的子域${\mathit{\boldsymbol{G}}_{\rm{1}}}$，其中绿色顶点为插值得到的内部新点；图 7(c)为更新区域$\mathit{\boldsymbol{G}}$，以及其边界与实顶点信息; 图 7(d)为确定第2次分解后得到的子域${\mathit{\boldsymbol{G}}_{\rm{2}}}$，并再次更新区域$\mathit{\boldsymbol{G}}$的边界与实顶点信息；图 7(e)为区域$\mathit{\boldsymbol{G}}$的边界实顶点中已无凹点，无须再分解；图 7(f)为最终子域分解结果。

本节将研究基于多目标优化的高质量网格生成方法，主要包括网格顶点间拓扑连接信息生成和网格顶点几何位置确定两个步骤。

由子域边界，生成每个子域内部的拓扑连接信息是本节所要解决的主要问题，也是生成高质量四边形网格的关键。本文主要采用文献[14]提出的拓扑信息构造方法，该方法可针对$N$边域生成内部网格拓扑信息(2≤$N$≤6)，即输入$N$边域中每条边所含有的折线单元数量，便可得出此种情形下所对应的顶点连接模板。该方法的主要思想是利用整数规划技术，使得所得到的拓扑连接模板具有最少的奇异点。因此，对于在第4节所得到的子域边界，首先要确定其$N$个角点的位置，然后根据简化对应原则计算其每一条边的顶点数量，再依据文献[11]中的模板生成其内部的顶点拓扑连接信息。

首先在子域边界中存在实顶点和虚顶点，对于实顶点即为该子域边界中的几何特征点。为进一步减少奇异点数目，采用4边域或3边域为基本区域单位。即为了能确定其角点信息，需要使得这些子域边界的实顶点个数保证为3或4。因而，将子域边界按照实顶点的个数$N$分为两类情形考虑：

1)$N>4$：对于具有4个以上实顶点的子域边界，需要对只由实顶点所构成的边界进行再次化简，使得边界实顶点个数为4。简化方法为第3节所述的边界化简方法。其化简参数可简化角度$A$设为90度，可简化面积比率$\gamma $设为0.25。并且免除原先的两个限制条件。这样简化到最后将只会保留到3个或4个顶点。

2)$N=3$或$N=4$：当边界只拥有3个或4个实顶点时，则无需进行化简，实顶点即为角点。

在角点确定之后，就可以计算其每一条边侧的边的数量，对于4边域，可以计算其4条边侧的边的数量$T, R, B, L$。对于3边域可以计算其3条边侧的数量$R, B, L$，根据文献[14]的方法把不同模式的参数计算完毕，就可以确定其边界内部的拓扑信息。对于$N=4$的情形，通过计算可知该区域需要填充的结构化网格的数量，但垂直填充和水平填充可分别在左右和上下两个边侧上任意填充，本文的实现原则是使拓扑模式中的非结构化网格位于整个区域的中间^[15]，即需要在垂直和水平两个相对方向填充数量基本一致的结构化网格。对于3边区域来说，$\alpha $的取值则为其变化范围内的平均值，使模式内的非结构化的网格大小适中。确定了这些参数之后，不同区域的内部网格拓扑信息便完全确定。对于$N=3$和$N=4$两种情形，文献[14]根据每条侧边的顶点数目，将其拓扑连接方式分别归结为2种模式和5种模式进行考虑。图 8中列出了一些模式的图例，对于其他的情形，可具体参考文献[14]中的拓扑连接信息的确定方法。

为生成高质量的四边形网格，利用多目标优化方法确定内部顶点的位置信息。首先，将文献[16]中的Laplacian算子推广到四边形网格，并将其度量引入到待优化的目标函数之中。四边形网格在顶点${v_i}$处的光滑性条件为

$ {\mathit{\boldsymbol{v}}_i} = \frac{1}{{\left| {N\left( {{\mathit{\boldsymbol{v}}_i}} \right)} \right|}}\sum\limits_{j \in \mathit{\boldsymbol{N}}\left( {{\mathit{\boldsymbol{v}}_i}} \right)} {{\mathit{\boldsymbol{v}}_j}} $

式中, $\mathit{\boldsymbol{N}}\left( {{\mathit{\boldsymbol{v}}_i}} \right)$为顶点${v_i}$的1-环邻域的顶点集合。由此，可定义四边形网格上的光滑性度量为

$ {F_{{\rm{smooth}}}} = \frac{1}{2}{\sum\limits_k {\left\| {\mathit{\boldsymbol{L}}{\mathit{\boldsymbol{v}}_k}} \right\|} ^2} $

式中，$\mathit{\boldsymbol{L}}$为Laplacian算子，其矩阵中的元素为

$ {L_{i, j}} = \left\{ \begin{array}{l} 1\;\;\;\;\;\;\;\;\;\;\;\;\;\;i = j\\ - \frac{1}{{N\left( {{\mathit{\boldsymbol{v}}_i}} \right)}}\;\;\;j \in N\left( {{\mathit{\boldsymbol{v}}_i}} \right)\\ 0\;\;\;\;\;\;\;\;\;\;\;\;\;其他 \end{array} \right. $

其次，为了得到面积大小尽可能均匀的四边形网格单元，从统计学及概率论中的方差概念可知，若所生成网格中的四边形单元面积的方差为零，则该网格即为单元大小均匀网格。因此首先计算该网格中所有4边单元的平均面积${A_{{\rm{ave}}}}$，即

$ {A_{{\rm{ave}}}} = \frac{1}{N}\sum\limits_{i = 1}^N {{A_i}} $

则平面网格上四边单元的均匀性度量可定义为

$ {F_{{\rm{uniform}}}} = \frac{1}{N}\sum\limits_{i = 1}^N {{{\left( {{A_{{\rm{ave}}}} - {A_i}} \right)}^2}} $

通常来说，理想的四边网格应满足在正则点处其4个夹角为90°，在奇异点处的$\rho $个夹角为$\cos \frac{{2\pi }}{\rho }$，其中$\rho$为该网格顶点的度数。因此，可统一定义相应的正交性度量为

${F_{{\rm{orth}}}} = \sum\limits_{i = 0}^N {\sum\limits_{k = 1}^\rho {{{\left( {\frac{{\mathit{\boldsymbol{T}}_k^i \cdot \mathit{\boldsymbol{T}}_{k + 1}^i}}{{\left\| {\mathit{\boldsymbol{T}}_k^i} \right\|\left\| {\mathit{\boldsymbol{T}}_{k + 1}^i} \right\|}} - \cos \frac{{2\pi }}{\rho }} \right)}^2}} } $

式中，${\mathit{\boldsymbol{T}}_k^i}$为在顶点${v_i}$处第$k$个边向量，$\mathit{\boldsymbol{T}}_{\rho + 1}^i = \mathit{\boldsymbol{T}}_1^i$。

综上，本文所研究的网格顶点几何位置生成问题可归结为一个多目标的非线性优化问题：即由3.4节中所构造的网格曲面${\mathit{\boldsymbol{M}}^{\rm{0}}}$作为初始解，优化变量为其内部网格顶点，优化函数为

$ \mathop {{\rm{arg}}\;{\rm{min}}}\limits_{{P_k}} \;{\omega _1}{F_{{\rm{smooth}}}} + {\omega _2}{F_{{\rm{uniform}}}} + {\omega _3}{F_{{\rm{orth}}}} $

式中，${\omega _1}, {\omega _2}, {\omega _3}$为加权系数。

对于上述优化问题，本文通过拟牛顿非线性优化算法进行求解^[17]。图 9给出了基于本文方法所构造的若干网格生成实例，并分别与文献[18]中的波前法，文献[7]中的铺砌法，以及文献[6]中的三角网格转换方法(见表 1)所生成的四边形网格生成结果进行了比较。在四边形网格模型中，顶点连接度不为4的内部顶点以及顶点连接度不为3的边界顶点称为奇异点。奇异点数目越少，表明四边形网格质量越高。由图 9可以发现，与文献[6-7, 18]中的3种经典方法相比，本文方法所构造的四边形网格的奇异点数目大大减少(见表 2)。

此外，为进一步验证本文方法的有效性，本文采用比例雅克比(Scaled Jacobian)作为四边形网格单元质量的评价标准^[19]。比例雅克比反映了网格的基本质量，如网格形状、内角大小等，其取值范围为[-1, 1]，越接近1则说明该网格单元质量越好，而比例雅克比的均方差越接近于零，则说明网格单元的大小越均匀。表 2列出了各个网格实例的相关数据，包括边界顶点数目，网格顶点数目，网格单元数目，奇异点数目，生成时间，比例雅克比的最小值，平均值和均方差。可以看出，在同样的离散边界下，本文方法与文献[6-7, 18]中的3种经典方法相比，所生成的四边形网格顶点和单元总数目较少，网格单元质量基本类似，计算时间成本大致相同，但奇异点的数目大大减少，更加满足工程仿真要求。

在计算效率方面，表 2中列出了各个方法相应的生成时间。在算法复杂性方面，本文方法与表 1中的3种经典四边形网格生成方法类似，时间复杂度均为O($N$)，其中$N$为给定的边界顶点数目。表 1对常用四边网格生成算法的主要特点和性能进行了总结与分析比较。

本文提出了一种基于边界简化与多目标优化的高质量四边形网格生成新框架。实例表明，与传统方法相比，本文方法生成的四边形网格具有较少奇异点，并具有良好光滑性、均匀性和正交性，可直接应用于有限元分析仿真求解。

将来的研究工作主要包括：

1) 由于该方法所生成的四边形网格无法从理论上保证没有翻转网格单元的存在，因而，如何对存在翻转单元的网格进行处理，是下一步工作需要研究的问题。

2) 将来可采用遗传算法求解本文中的多目标优化问题，进一步优化网格质量。

3) 进一步将本文提出的方法应用于等几何分析中复杂计算域的参数化问题^[9]。