正交六面体网格生成算法

我爱前处理 2023年3月7日浏览：1501 收藏：2

1. 问题描述

在用有限元法或者有限体积法求解流体力学问题时，需要先将求解区域划分成网格。区别于在物体表面生成的网格（surface mesh），我们称这种划分三维区域的网格为体网格（volume mesh）。

体网格根据其单元形状可以分为四面体网格（tetra-mesh），六面体网格（hexa-mesh），以及四面体或六面体为主的多面体网格（tetra/hexa-dominated mesh）。而根据其生成方式又可以分为结构化与非结构化网格（stuctured, non-structured），贴体与非贴体网格（conformal, non-conformal），等等。

对于四面体网格，较为常用的生成算法包括Delaunay法和波前法（advancing front）等。而对于六面体网格，较为常用的算法有：映射法（mapping），扫掠法（sweeping），以及正交切割单元法（Cartesian cut-cell）。映射法和扫掠法只适用于特定类型的几何模型；而正交切割单元法具有较强的普适性，只需要提供模型的表面网格就可以自动生成六面体为主的多边形网格（并非纯六面体网格）。这里我们主要介绍一下这种方法。

2. 求解流程

有些体网格生成算法直接从描述几何模型的参数方程出发（例如映射法，扫掠法），而另一类体网格生成算法从模型的表面网格出发（例如Delaunay法，波前法），正交切割单元法属于后者。它对于输入的面网格有一些要求：

1) 纯三角（triangular）：所有单元均为三角形。

2) 满足水密性条件（water-tight）：任意一条边恰好属于两个三角形单元。

3) 满足流形条件（manifold）：任意两个三角形要么不相交，要么共享一条边。

算法最终得到的是一组网格单元（cut-cell element），每个网格单元由一个正方体和若干个切割面（cut-face）共同构成。在必要时可以计算出正方体被切割后形成的多面体，所有这些多面体共同构成模型内部空间区域的一个剖分。

算法的大致流程如下：

1) 初始化：首先根据面网格的包围盒（bounding box），以及用户指定的分辨率（resolution）生成均匀的正方体单元，作为第一层网格单元的基础。

2) 切割：对于输入的面网格中的每个三角形，确定所有与其相交的正方体单元，并计算出被正方体切割后得到的切割面（cut-face），储存在相应的正方体单元中。每个三角形上的所有切割面构成这个三角形的一个剖分。

3) 细分：对于在模型内部而未与表面网格相交的正方体单元，可直接放入体网格当中。对于与表面网格相交的正方体单元，根据其存储的切割面数量和法向，决定是否将其细分（refine）。对于满足细分条件的单元，将其分成大小相等的8个小正方体单元，并将原本单元中存储的所有切割面再次切分后，分别存储到对应的小单元之中。

4) 组装：当所有正方体单元均不满足细分条件，或者达到预设的最大细分次数时，细分过程终止。将所有含有切割面的单元放入体网格当中。此时，体网格中包含了所有完全位于模型内部的单元，以及与模型表面相交的单元，它们构成模型内部空间区域的一个剖分。

3. 数据结构

为了完整存储一套体网格的全部数据，我们需要记录三个列表：节点坐标的列表（point list），切割面的列表（face list），以及网格的单元列表（cell list）。这里我们把节点和切割面单独存成列表，而不是放到网格单元里面，是为了避免重复记录相同的数据。例如一个切割面属于两个网格单元，那么我们只需要记录一份这个面的数据，然后在两个网格单元中分别存储这个面在列表中的下标即可（相当于记录了它的地址）。

切割面（cut-face）的存储切割面是由若干个节点构成的平面多边形，它是从输入面网格中的某个三角形被正方体切割得到的。对于每个切割面，我们只需要记录构成它的所有节点，以及它所在平面的法向量即可。切割面（cut-face）的存储切割面是由若干个节点构成的平面多边形，它是从输入面网格中的某个三角形被正方体切割得到的。对于每个切割面，我们只需要记录构成它的所有节点，以及它所在平面的法向量即可。

网格单元（cut-cell）的存储

每个网格单元由一个正方体和若干切割面构成，在需要时可以将其转化为多面体单元。由于自动细分机制，一个单元若被细分，则存储其 8 个子单元的地址（下标），因此全部网格单元构成一个八叉树结构。

如图，这个网格单元中包含 4 个切割面，每个切割面由面网格中的一个三角形和正方体相交得到。它们将正方体切割成两个多面体。只要记录这些切割面和立方体，就可以在需要的时候计算出切割后形成的两个多面体。

4. 计算方法

这里我们重点讨论第一步（切割）中用到的计算方法。细分和组装的部分我们这里先不详细讨论（如有需要后面可以为此另写一篇专题）。切割三角形的计算可分为三步：

1）对于每条三角形边，计算其与正方体平面的交点，放入节点列表中；同时，将节点有序插入到这条边的节点列表中。这里“有序”指的是每个节点的前，后节点为几何上与其相邻的节点。在 C++中可以借助 std::map 来实现有序插入。

2）将第一步算出的交点两两组成混合边（hybrid-edge）。对于每条混合边，计算其与正方体平面的交点，放入节点列表中；同时，将节点有序插入到这条边的节点列表中（这部分与第一步类似）。

3）将第一步和第三步生成的节点组装成平面多边形（即切割面），放入切割面列表中；同时，将切割面的下标记录在对应的网格单元中。

下面具体讨论这些步骤中用到的特殊方法和结构。

关于混合边：

混合边指的是三角形与切割它的平面的交集所在的线段。由于三角形与平面均是三维空间中的凸集（convex set），因此它们的交集也是凸集，这意味着无论它们相对位置如何，交集都只有一条线段（包括退化为一个点或者是空集的情况）。也就是说，只要我们记录下每个与三角形相交的平面上产生的交点，最后就能得到一条在这个平面上的混合边。在切割三角形的计算中，混合边与三角形边有着同等的重要性（统称为“边”）。

有两种退化情况需要注意：第一种是如果三角形的一条边恰好落在某个平面上，那么可以跳过计算与平面相交的过程，而直接把这条三角形边作为平面上的混合边；另一种情况是，如果某个平面与三角形的边相交得到的两个交点过于接近（小于浮点数误差），或者恰好经过三角形的某个端点，那么应当将这个端点作为这个平面上退化形式的混合边。

另外还需要注意的是，在计算混合边与正方体平面的交点时，若不加处理，会导致重复计算（见下右图）。因此对于这种交点，需要记录其所在的网格线（grid line），确保每条网格线上只有一个这样的节点。如果计算新的交点时发现对应的网格线上已经存在一个这样的节点，那么应将这个节点作为交点，而不是生成新的节点。

关于组装多边形：

在将节点组装成平面多边形的过程中，我们也需要用到一些计算方法。

1）首先，我们需要记录每个节点的相邻节点（neigh bor），这可以通过按顺序遍历每条边（三角边或者混合边）上的节点来得到。

2）接下来，对于每个节点，又需要将其所有相邻节点按照逆时针方向排序（因为它们都在同一个平面上，所以是可以做到的）。关于这部分的具体算法在参考文献[3]中有详细的描述，这里就不再赘述了。

3）我们将所有未被组装的节点有序对（即半边，half-edge）放在一个集合当中。从一条半边出发，根据相邻节点的顺序找到下一个节点，直到回到初始节点，得到一个完整的（顺时针）回路构成一个平面多边形。将这个过程中访问到的半边从集合中移除，直到集合为空，此时我们便得到了三角形上所有的平面多边形。

5. 一些结果

这里展示一些用我们自己的代码生成的正交网格结果。我们将结果转化为多面体网格的形式（即计算出每个网格单元中切割后的多面体，并将位于模型内部的多面体放入一个列表中），并将其输出为 Paraview 可以打开的文件格式。为了更直观的看到模型内部单元的形状和分布，我用 Paraview 中的 clip filter 工具切掉了网格中的一些多面体单元。由于还未做细分的部分，这里显示的仅为用均匀背景网格切割得到的结果。