四边形网格的形态好，计算精度高。四边形网格的生成严重依赖于几何形状，对简单的平面区域生成较为简单；对复杂的区域，，可以先将区域分割为简单的图形，再划分。有如下方法：有限元差值法，扫略法，合并法。

1.插值方法（FE Interpolation）

插值方法是利用插值思想，计算出区域内，要离散四边形的的顶点坐标。插值思想与有限元四边形单元的插值函数非常类似，参考《形函数》。对于下图所示的，任意平面内的四边形（四条边都是直线）做FE插值。先将任意四边形推广当广义坐标系下[-1,1]，以给定尺寸，计算每个方向上的四边形。

推导公式如下：

编码实现的部分代码如下：

FEInterpolation::FEInterpolation( Vertex& v0,Vertex& v1, Vertex& v2, Vertex& v3, double len) { //以四条边中，相邻的两条最长边作为划分个数依据 double dis01 = v0.distance(v1); double dis12 = v1.distance(v2); double dis23 = v2.distance(v3); double dis30 = v3.distance(v0);
double maxU = std::max(dis01, dis23); double maxV = std::max(dis12, dis30);
int u = int(maxU / len); int v = int(maxV / len);
double uParaInterv = 2.0 / u; double vParaInterv = 2.0 / v;
for (int i = 0; i < u; i++) { double epsi = -1 + i * uParaInterv; for (int j = 0; j < v;) { //计算H值 double eta = -1 + j * vParaInterv; double H1 = (1 - epsi)*(1 - eta) / 4; double H2 = (1 + epsi)*(1 - eta) / 4; double H3 = (1 + epsi)*(1 + eta) / 4; double H4 = (1 - epsi)*(1 + eta) / 4;
//v0 Vertex m_v0 = v0 * H1 +v1 * H2 + v2 * H3 + v3 * H4;
//计算v1, i+1, i++; epsi = -1 + i * uParaInterv; H1 = (1 - epsi)*(1 - eta) / 4; H2 = (1 + epsi)*(1 - eta) / 4; H3 = (1 + epsi)*(1 + eta) / 4; H4 = (1 - epsi)*(1 + eta) / 4; Vertex m_v1 = v0 * H1 + v1 * H2 + v2 * H3 + v3 * H4;
//计算v2,j+1 j++; eta = -1 + j * vParaInterv; H1 = (1 - epsi)*(1 - eta) / 4; H2 = (1 + epsi)*(1 - eta) / 4; H3 = (1 + epsi)*(1 + eta) / 4; H4 = (1 - epsi)*(1 + eta) / 4; Vertex m_v2 = v0 * H1 + v1 * H2 + v2 * H3 + v3 * H4;
//计算v4,i--; i--; epsi = -1 + i * uParaInterv; H1 = (1 - epsi)*(1 - eta) / 4; H2 = (1 + epsi)*(1 - eta) / 4; H3 = (1 + epsi)*(1 + eta) / 4; H4 = (1 - epsi)*(1 + eta) / 4;      Vertex m_v3 = v0 * H1 + v1 * H2 + v2 * H3 + v3 * H4; QUAD quad(m_v0, m_v1,m_v2, m_v3);      quads.push_back(quad); }  }}

下图是划分的结果：

2. 映射（mapping）

上述方法只能处理边界是直线的情况，对于边界是曲线的方法，需要通过扫略方法。对于如下图所示的，四条边界为曲线的四边形，建立映射方程。可以发现，此映射方程同时考虑角点和曲线弯曲情况。

此方法与章节1的实现类似，以下是通过这种方法划分的网格。

对于某些简单结构，可以采用扫略（Sweeping）方法。其核心还是计算四边形节点坐标，只是不再采用插值函数计算；而是，以初始点，沿着切向量的增量计算新的点。如下，左图的情况较为简单，直接沿着路径方向扫略生成新的四边形；右图所示的两条扫略路径为曲线，沿起始点的切线（任意一条切线的路径均可）方向扫略，形成新的点后，将两段的点向曲线投影，并对新的点做距离规整处理；

3. 合并三角形为四边形

对于复杂区域，直接划分四边形难度很大。一般先划分为非结构三角形，然后将三角形合并为四边形网格。对于三角形，要能够快速访问其相邻三角形的信息，参考《三角形数据结构》。有如下技术问题要解决：

• 三角形与那个三角形合并
• 那些三角形合并，而其它的三角形不合并

这两种都可以通过衡量三角形与四边形质量来解决。

3.1 质量参数

三角形的形状通过α参数（α-quality coefficient）定义，如下式：

四边形的扭曲系数，通过下式定义：

典型四边形的扭曲系数如下图：

临界值γ定义为允许合并的最小值。

3.2 合并算法

合并算法的步骤如下：

• 计算三角形每条边的β值，如果某条边没有相邻三角形，β值设定为-1，意味着此边不能删除，是边界。

• 合并，依次查看三角形的每条边，如果β值大于临界值γ，合并。并且将合并前三角形的其他四条边的β值赋值为-1（因为四边形的边界再不能被删除）；

• 继续合并，直到所有边都不能被合并。
  

此算法的主要任务是计算每条边的β值，其余步骤都是很清晰的。具体实现后续再提供。

3.3 临界值γ的影响

γ值越小，合并的三角形就越多；如下图所示，

4. 直接法

4.1 区域分割

将区域分割为，能满足第1,2节所示的几何，划分好子几何后，拼接起来即可。

4.2 铺层方法

此类方法是采用类似AFT的思想，以一条边界为前沿，搜索四边形。难点在于边界的填充处理。

4.3 背景网格方法

采用四叉树管理区域，此时不再是三角形，而是直接生成四边形。难点在于边界、过渡处网格的处理。

当然，还有基于标架场和黎曼度量的方法。

路漫漫其修远兮，共勉 。

参考：

Finite Element Mesh Generation