三角形网格算法

① 三角曲面求交的网格法

6.1.3.1 网格法的基本原理

由于三角曲面求交运算最终是三角形之间的求交运算，因此，我们在不能漏掉任何可能相交的三角形对的同时，要尽量排除不可能相交的三角形对。网格法的目的就是寻找可能相交的三角形对。该方法先将曲面包围盒的相交区域划分成规则六面体格网，先找到与每个网格单元相交的三角形，再将每个单元内分别来自两张曲面的三角形组成可能相交三角形对，再对可能相交三角形对进行求交运算。网格法的基本过程为:

(1)先在合适的坐标系下分别构造出两张三角曲面的矩形包围盒，求出两个包围盒的交作为新的包围盒;

(2)然后将新包围盒划分成由若干规则六面体组成的格网;

(3)针对每个网格，分别从两张曲面中找到与其相交的所有三角形，并组成可能相交的三角形对;

(4)对所有网格中的三角形对进行整理，删除重复三角形对;

(5)对所有三角形对进行求交运算，计算出交线段;

(6)进行交线追踪，获得曲面交线。

6.1.3.2 格网的构建

为了方便，常常将两张曲面的包围盒的相交区域划分成等尺寸的规则六面体格网。这样，仅仅根据坐标就可以很方便地确定每个结点位于哪个网格内。

显然，不与任何一个相同网格相交的两个三角形是不可能相交的。网格法可以排除很多不相交的三角形对，也不会漏掉任何相交的三角形对。网格的大小是影响计算效率的关键。如果网格单元尺度太大，则会得到更多的可能相交三角形对，降低了不相交三角形对的排除效率;如果网格单元尺度太小，内存就会开销增大，而且三角形与格网相交计算量也会增大。一般情况下，可以使网格边长与三角网格的平均边长相当。

6.1.3.3 三角形与格网的相交判断

由于曲面上有数量众多的三角形，如果直接判断每个三角形与每个网格单元是否相交，则计算效率非常低。下面介绍一种判断一个三角形与格网相交的方法，算法思路如下:

(1)在格网所在坐标系下，将三角形与格网分别投影到xoy面上(图6.3(a))。格网的投影为一个平面格网，而三角形的投影为三角形或一条线段。如果三角形的投影为线段，则选择xoz或yoz为投影面。为了叙述简便，这里只介绍在xoy面上的投影为三角形的情况。

(2)计算投影三角形的边与平面网格线的交点，并计算过每个交点的xoy面的垂线与三角形的交点，将这些交点与三角形的顶点合并组成点集Ω。

(3)利用Ω中的点在xoy面上的投影坐标直接判断每个点所在的格网单元(包括位于格网单元内部或边界上)，找到同属一个格网单元内的点组成子集(一个点可以位于多个格网单元内)。

(4)针对每个Ω的子集Ф及投影所在的格网单元g(图6.4(b))，将原三维格网中平面投影为g的所有格网单元组成列G。找到子集Φ中点的最小与最大z坐标，从最小到最大z坐标所跨越的G中的网格单元即为与该三角形有相交的格网单元。

(5)删除由上述方法所得的重复记录的格网单元，就可得到与该三角形相交的所有格网单元。

图6.4 三角形与格网求交示意

6.1.3.4 网格法曲面求交算法

网格法三角曲面求交算法如下:

三维地质建模方法及程序实现

② 推进波前法

经过近年来的发展，推进波前法(Advancing Front Technique，AFT)已经成为通用的全自动非结构化有限元网格生成方法之一。该方法最初由Lo提出并用于平面区域三角形网格自动生成(Lo，1985)。AFT方法具有生成的边界网格质量高、易于自适应加密等优点，但不同于Delaunay三角剖分算法，它没有后者那样成熟的理论依据，在很多情形下靠经验解决问题，但是这并不妨碍它的成功应用。

推进波前法的基本思路是:按照剖分规模将边界离散成有序线段，然后从边界出发，依次以边界线段为三角形的一条边，在边界点与内部点中寻找合适点，组成三角形，选取组成三角形顶角最大的点为最终三角形顶点;将已形成三角形的边界线段从边界链表中删除，形成新边界;重复上述过程直到除边界外的三角形的边两侧均有三角形为止。为了更好地说明该算法，下面先介绍几个术语。

3.3.1.1 二维AFT方法术语定义

(1)剖分域:即需要剖分的区域。正确地定义剖分域(区域的几何描述)是网格能够正常生成的必要条件。剖分域是由一系列有向边界曲线围成的连通域，并且每条边界曲线必须是简单封闭曲线。通常情况下，剖分域的外边界按照逆时针排列，而内边界则按照顺时针排列。

(2)前沿Ω:所有未剖分区域的边界线段以及端点的集合构成Ω，前沿包括活跃前沿(记为Ω₁)和非活跃前沿(记为Ω₂)，其中活跃前沿为当前正在推进的前沿，非活跃前沿为暂时不推进的前沿。

(3)选定前沿S:选定前沿S是Ω中的一个元素。S的选取取决于网格的生成策略，如果为了保证生成网格的尺寸过渡以及保证小尺寸单元优先生成，一般选取Ω中的最小前沿作为S。如果为了程序实现上的便利，则从Ω中从前往后依次选定一个元素作为选定前沿S。

3.3.1.2 算法要点

(1)选取合适的数据结构建立点、边、三角形之间的关系，并建立储存点、边与三角形的链表。

(2)选取合适的驱动方式。如以三角形的边为基础进行波阵式扩展，必须考虑边的使用次数与方向:任何位于区域边界上的边应且只应使用1次，任何位于区域内部的边应且只应使用2次(正、反方向各1次)。因此，在初始状态，应将边界边的使用次数赋1，内部边使用次数赋0。

(3)以边为基础进行波阵式扩展，是以某边为三角形的一条边，再从点集中寻找合适的顶点组成三角形的过程。所寻找到的点必须满足以下要求:新形成边与已生成的边不能相交;所有边必须满足使用次数要求(边界边使用一次，其余边使用两次);新顶点与该边(有方向)组成的三角形面积必须大于零;保证顶角最大。

3.3.1.3 算法与程序代码

平面区域的AFT方法主要有三大步:向剖分域中布点、离散剖分域的边界和推进前沿生成三角形。

3.3.1.3.1 布点

布点即是根据需要得到的三角形单元的各边的大概长度，在剖分域内生成一系列的散乱点。最常用的方法是先根据剖分域边界上端点的x和y坐标的最大值和最小值，生成一个包含剖分域的矩形，该矩形也叫做剖分域的包围盒;再在矩形中生成点，最简单的是生成“棋盘状”的一系列点，另外是生成“正三角形状”的一系列点;生成一系列点之后，判断这些点是否落在区域内，若是，则为需要布设的数据点，否则，删除;同时需要注意的是若某些点落在了区域内，但是又距离边界太近，依旧删除这些点。图3.9中(a)为剖分域，(b)则是布设“正三角形状”点的结果。

图3.9 在平面区域中布点

3.3.1.3.2 离散边界

离散边界即是按照剖分规模或需要得到的三角形单元的各边的大概长度将边界离散成有序线段，如图3.10所示，为布点之后进行离散边界的结果。

图3.10 离散边界

3.3.1.3.3 生成三角形

以三角形的边为基础进行波阵式扩展生成三角形，即以三角形的边为推进前沿，主要过程有如下4小步。

第一步:建立点集PS和边集ES。初始点集PS包括所有布设的数据点和边界离散后小线段的端点。初始边集ES只包括边界离散后的有向线段。此时，边集ES就是前沿Ω。

第二步:以边集ES中的边E_i为基础搜索顶点P_i，即选定E_i作为选定前沿S，以该点为顶点、该边为一边形成三角形。设E_i的端点为A与B，所有待选点P_i与E_i组成的顶角为∠AP_iB，将顶角从大到小排序，从最大顶角开始，依次选择对应的顶点P_i与E_i组成三角形。如果形成的三角形满足以下要求，则为新三角形，P_i为合适的顶点:①新形成三角形的边与已生成的边不能相交;②所有边必须满足使用次数要求;③新顶点与该边(有方向)组成的三角形面积必须大于零。如果不满足则选下一个顶点。

第三步:找到合适顶点后，将新顶点与选定前沿S(即边E_i)的端点连成的边加入到边集中，生成新的前沿Ω的元素，并将新形成的三角形加入到三角形集中，删除原选定前沿S，选定边集ES中的下一条边作为选定前沿S。

第四步:重复第二、三步，当所有边满足使用次数要求时循环结束。图3.11中，(a)、(b)和(c)依次为循环一步、二步和多步之后形成的三角形，当循环结束时得到的三角剖分如图3.12所示。

图3.11 推进前沿过程

三维地质建模方法及程序实现

函数CreateTrgls()中调用的函数SearchID()搜索某一点在顶点集合surf-＞pNodes的ID;函数CountCos()用于计算当一条线段/边搜索到一顶点并组成三角形时该顶点处角度的余弦值;函数SameLine()用于判断两条线段/边是否相同。

3.3.1.4 约束的处理

区域三角剖分中的约束是指待剖分区域中存在特定的点或线，称为点约束或线约束，其中约束点必须是剖分后网格的顶点，而约束线必须是剖分后网格三角形的边的集合，不存在某个三角形跨越约束线的状况。

3.3.1.4.1 点约束的处理方法

点约束的处理非常简单，直接将约束点加入到生成的点集中，再删除与约束点距离非常近的点，然后就可以按无约束的方法进行三角剖分了。

3.3.1.4.2 线约束的处理方法

当待剖分区域存在线约束时，可以将所有线约束作为一种边界。在剖分前，与外边界及内边界的处理方法一样，先按照一定的规模将约束线离散成顺序连接的线段，每条线段均作为三角网格的边，然后设置约束线上的边的使用次数为0，并加入到原始边集中，再按照按无约束的方法进行三角剖分即可完成约束三角剖分。图3.13(a)为含线约束的待剖分区域，图3.13(b)为约束剖分网格。

图3.13 约束三角剖分实例