有限元网格是有限元走向工程应用的桥梁，在有限元分析中占有重要的地位。有限元网格的生成是指将连续的几何模型离散成为计算所需的单元信息，其质量的优劣直接影响求解的速度与精度。本文简要为大家介绍几种常用的有限元网格生成方法。
01 映射法 (Mapping)
映射法是最早出现且最直观的一种有限元网格生成方法。在有限元网格发展早期，常采用扫掠法和拖拽法对简单四边界区域进行网格生成。将已剖分的起始边上各点沿特定的方向移动若干次至终止边，即完成对区域的剖分。此类方法实际上就是一种映射，采用向量作为坐标变换方式。
扫掠法
拖拽法
映射法的核心在于坐标变换，其基本步骤如下：
- 通过适当的映射函数将待剖分物理域映射到参数空间中形成规则参数域；
- 对规则参数域进行网格剖分；
- 将参数域的网格反向映射回物理空间，从而得到物理域的网格。
映射函数有很多种，可根据实际问题进行选取。若区域的四条边都为直线，可以采用四节点四边形的拉格朗日插值函数进行坐标变换；若区域的边界为曲线，则可以将边界表示为参数形式C=C(t)，同样可以将待剖分物理域映射至单位正方形参数域。此外，还可以采用偏微分方程法，通过数值求解偏微分方程而得到参数空间与物理空间之间的映射关系。
拉格朗日插值函数
将边界表示为参数形式C=C(t)
映射法最大的优点在于既能生成结构化网格，也能生成非结构化网格。映射法对单连通区域算法简单、速度快、单元质量好，并且适用于曲面网格生成，可与形状优化算法集成。因此，映射法在众多的商业有限元分析软件中占有重要的地位。
然而，映射法在实际应用中也存在许多难点，首先便是如何将复杂的多连通区域分解为若干可映射的子区域。目前在实际应用中，采用映射法处理复杂几何体的时间往往远大于网格生成的时间。此外，处理子区域之间网格的协调性和指定区域网格的疏密过渡问题，对用户有着较高的要求。
02 四（八）叉树法 (Quadtree/Octree)
四（八）叉树法是一种基于栅格法的全自动非结构化网格生成方法。栅格法的特点是用一组背景栅格覆盖剖分区域，删除区域之外的栅格，并以剩余部分作为初始条件生成网格。Yerry和Shephard首先将用于近似表达几何对象的四（八）叉树空间分解法引入到网格剖分领域，利用四（八）叉树数据结构，可对背景栅格递归细分，逼近边界。
以二维为例，四叉树法的基本步骤如下：
- 将剖分区域的边界离散化，得到一组线段集和点集；
- 创建背景栅格覆盖目标区域，根据点集的分布平衡四叉树；
- 保留与目标区域相交的栅格，删除完全落在目标区域之外的栅格；
- 内部栅格与边界栅格的相容网格剖分；
- 网格优化。
四叉树法
四（八）叉树法的难点在于内部栅格与边界栅格的相容网格剖分。以四叉树为例，在平衡叉树的步骤中，通常会将相邻区域之间的细分等级之差限制在1以下，使每个区域仅会在顶点和边中点处产生网格点，以减少相容网格剖分时所需的模板数量。然而，对于三维八叉树而言，即使相邻区域之间细分等级之差不大于1，也需要建立至少78种模板。为此，有许多学者致力于降低八叉树法程序实现的复杂度。
四叉树相容剖分模板
03 Delaunay三角化法(Delaunay Triangulation)
谈到Delaunay三角化的来历，就不得不先提它的对偶图——Voronoi图。1908年，俄国数学家Voronoi提出：给定平面内点集Q，对于其中每个点p，存在一个区域，该区域内的任意一点与p的距离是与Q中所有点距离的最小值，这些区域的合集构成Q的Voronoi图。1934年Delaunay提出，通过连接Voronoi图相邻区域的点，可以得到一个平面点集的三角化，这样的三角化就称作该点集的Delaunay Triangulation，简称DT。
Voronoi图与Delaunay三角化
给定一个平面点集Q，其DT有以下特性：
- 唯一性：在不考虑四点共圆的情况下，DT是唯一的；
- 最大化最小角：在Q的所有三角化中，DT使最小角达到最大；
- 空圆特性：DT中每个三角形单元的外接圆都不包含其它节点，此特性也称为DT条件，可用于判断三角化是否为DT，如下图所示。此外，通过简单的边翻转（BD->AC）可将局部非DT转化为DT。
边翻转示意图
也就是说，DT避免了生成小内角的长薄单元，尽可能地将所有三角形趋向均匀。鉴于其良好的特性，到了上世纪80年代，DT开始大量应用于有限元的网格生成，出现了许多DT的生成算法，最常见的当属Bowyer-Watson算法。Bowyer-Watson算法采用逐点插入法，每次插入新点后，都要保证当前的三角化符合DT条件。此外，点集的插入顺序是影响生成算法效率的重要因素，一般采用随机次序插入。Bowyer-Watson算法的主要步骤如下：
- 插入一个新点，在已生成的三角形中找出所有外接圆包含插入点的三角形，得到影响区域多边形；
- 删除多边形的内部边，连接插入点与多边形的全部顶点，连接前需要进行可见性检查。
Bowyer-Watson算法示意图
DT生成算法对凸区域十分有效，然而待剖分区域很有可能在局部是凹的，凹区域采用DT方法生成的网格不能保证边界的完整性，需要进行边界恢复。对二维而言，可以通过简单的边翻转操作恢复边界。若AB是丢失的边界，循环多边形的每对相邻三角形构成的四边形进行边翻转判断，确保四边形的对角线不会穿过AB，直至多边形内的所有三角形都不穿过AB即可。
二维边界恢复示意图
综上，二维的Delaunay三角化法的基本步骤如下：
- 得到待剖分区域的点集，包括边界节点和内部节点；
- 将点集中的点生成一个随机插入次序；
- 创建一个能包含点集凸包的三角形或正方形；
- Bowyer-Watson算法，确保每一步都满足DT条件；
- 通过拓扑操作恢复边界。
Delaunay三角化在由二维推广至三维的时候却遇到了一些问题。首先，三维Delaunay三角化会产生四个节点几乎共面的薄元，薄元的单元质量很差却能符合三维的DT条件。此外，三维的边界恢复不能直接由二维推广，需要分别针对边界边和边界面进行处理。不过随着研究的深入，这些问题目前都已经有了较好的解决办法。
04 前沿推进法 (Advancing Front Technique)
前沿推进法也叫波前法，由Lo和Lohner提出，是一种基于边界的全自动非结构化网格生成方法。顾名思义，该方法将区域边界作为初始波前，根据一定规则向区域内部推进，推进的过程中进行网格生成，当所有波前在区域内部聚合消失，网格生成完毕。与Delaunay三角化法具有成熟的理论依据不同，前沿推进法的推进规则在很多情形下是依靠经验制定的，但这并不妨碍它的成功应用。前沿推进法的算法时间复杂度与四（八）叉树法和Delaunay三角化法相当，并且生成的网格在边界处有更好的质量。
以二维为例，前沿推进法的基本步骤如下：
- 将待剖分区域边界离散化，得到一组线段集合，即初始前沿；
- 根据前沿各线段的位置关系和推进规则生成新节点，并与前沿线段连接得到新的单元，并更新前沿中的线段；
- 重复步骤2，直至前沿为空，剖分完成。
前沿推进法示意图
前沿推进法的优势在于新节点和新单元同时生成，便于利用已存在的单元控制后续生成的节点位置，以达到质量控制、局部加密或网格过渡的要求。此外，从边界向内推进，能够优先保证边界附近的单元质量。
不过，前沿数据结构的设计一直是此方法的难点。一方面要考虑推进规则的完整性；另一方面，在生成新单元时要进行大量的相交判断、距离判断，需要尽可能减少执行判断的次数以提升剖分效率。
二维前沿推进法的一种推进规则
就算法的时间复杂度而言，四（八）叉树、Delaunay三角化法和前沿推进法这三种全自动的非结构化网格生成方法计算效率都为O(nlogn)，难分伯仲。但实际上各种网格生成方法都有各自的优势和局限性，在实际的网格生成算法中，往往需要将多种方法组合运用，扬长避短，以达到更好的剖分效果。
图片出处
[1] S.H. Lo, Finite Element Mesh Generation, CRC PressTaylor & Francis USA, 2015
[2] PJ Frey, L Marechal, Fast Adaptive Quadtree MeshGeneration, 7th International Meshing Roundtable, 1998
