在计算流体动力学（CFD）中，空间离散化——具体而言即网格生成与分区过程——对数值模拟的精度和效率起着至关重要的作用。网格定义了连续物理域如何被划分为离散控制体，并在其上求解控制方程。总体而言，工程应用中主要采用两类网格：结构化网格与非结构化网格。
结构化网格具有规则的拓扑结构，其中每个内部节点拥有数量和排列一致的相邻单元，通常构成六面体单元。这种均匀结构允许基于网格索引隐式定义节点连接关系，无需显式存储邻居信息。结构化网格由于单元质量较高且网格间距变化平滑，通常具有更高的数值精度和更好的收敛性。然而，针对复杂几何体生成结构化网格往往费时费力，需要大量人工干预，且在许多情况下需对原始几何体进行简化，以实现拓扑规则的网格。
相比之下，非结构化网格缺乏这种规则的连接模式。内部节点可能拥有数量不等的相邻单元，这些单元可以呈现多种形状，如四面体、棱柱、锥体或一般多面体。非结构化网格的灵活性使其能够对几何复杂的域实现自动网格生成，适用于涉及复杂边界或多个相互作用组件的仿真。然而，这种灵活性也带来了内存需求增加和计算开销增大的代价，因为在数值计算过程中必须显式存储并访问邻居关系。
八叉树网格
八叉树作为一种组织空间对象的方法，基于单元密度对空间进行划分，从而避免在网格生成过程中遍历所有对象。其原理相对简单：当满足划分条件时，将三维空间等分为八个子空间，并将空间对象分配至相应子空间。八叉树网格作为一种常见的非结构化网格，其生成方式为：首先将计算域划分为一个或多个较大的立方体网格，再对这些立方体网格进行递归细分，分裂为八个子网格，直至每个子网格满足预定的尺寸要求或被几何边界裁剪为止。
尽管八叉树属于非结构化网格，但其内部不同尺寸的六面体单元之间保持固定的细化关系。这使得在存储网格时，可以间接存储网格单元的几何信息和邻居关系。本文介绍一种利用空间填充曲线（SFC）存储网格单元的方法。
空间填充曲线（SFC）
空间填充曲线是一个一维区间，包含从 n 维空间到一维空间的映射。从概念上看，它可以视为一条以特定顺序访问空间域内所有离散点的连续路径。这种顺序本身具有一定的空间位置特性，从而有助于优化缓存布局。
Morton 曲线是一种常见的空间填充曲线，将 n 维空间映射为有序的线性序列。应用于坐标时，Morton 编码定义了一条 Z 形空间填充曲线，因此也称为 Z 曲线。该编码确保在物理空间中邻近的点通常具有相近的 Morton 编码，这一特性对内存布局和并行数据访问大有裨益。
2.1 Morton 编码的生成
Morton 编码的生成本质上是对各维度坐标进行位交织操作。
以二维平面中最大细化层数 maxLevel = 2 的四叉树网格为例进行说明。树的根节点代表网格的包围盒，叶节点对应各个网格单元。每个单元基于其质心位置被分配一个六位二进制编码。Morton 编码的总位数为：
其中 numDim 表示网格的空间维度，maxLevel 表示最大细化层数。
每个四叉树单元的位编码可由其在笛卡尔坐标系中的相对索引确定。例如，图 3 中红色单元的 Morton 编码可表示为 (10,01,01)。代表初始最粗网格的网格边界 00 不纳入 Morton 编码。二进制 Morton 编码在需要时可转换为十进制；例如，红色单元的十进制序号为 37。该四叉树网格的所有单元可排列成如图 4 所示的树形结构，其中彩色单元对应图 3 中的二维四叉树单元。
生成所有单元的 Morton 编码后，从左下角单元遍历至右上角单元，可得到图 5 中黑色曲线所示的空间遍历序列。该曲线呈现出明显的锯齿形结构，因此 Morton 曲线也称为 Z 曲线。可以观察到，四叉树中的兄弟单元在 Morton 曲线上同样占据相邻位置。此外，在四叉树动态细化过程中，新生成的兄弟单元聚集在其父单元在 Morton 曲线上位置的附近。在三维空间中对八叉树进行细化时，八叉树单元的 Morton 编码也展现出类似的空间特性。
2.2 细化层数
生成 Morton 编码后，每个单元的细化层数可直接从 Morton 编码中推导得出。
对于维度为 dim 的 Morton 编码，算法如下：给定一个长度为 N 的 Morton 编码数组 mortonKey，以及一个八叉树/四叉树网格（即 Morton 编码的集合）。
例如，图 3 中粉色单元的 Morton 编码为 (11,00,00)。初始设 level = N/dim = 3。从 Morton 编码的第一层开始，编码 11（非零）跳至下一层。第二层的 Morton 编码为 00，将 00 的 x 位取反得到 10。检查 (11,10,00) 发现其已存在于四叉树中，表明在第 2 层存在兄弟邻居。继续至第三层，Morton 编码为 00，将 00 的 x 位取反得到 10，检查发现 (11,00,10) 不存在于四叉树中。因此 level = level - 1 = 2。由此可知，该单元具有两个细化层数。
2.3 单元几何信息
由 Morton 编码以及根节点单元的几何中心和边长，可以得到细化单元的边长、几何中心和顶点位置。
2.4 单元邻居信息
单元邻居信息对于偏微分方程（PDE）数值计算的重要性不言而喻。邻居信息的存储与检索对 PDE 求解器的整体计算效率有重大影响。对于采用 Morton 编码的四叉树/八叉树单元，邻居信息天然嵌入在 Morton 编码中，通过简单的位运算即可获取。
四叉树/八叉树单元的邻居分为兄弟邻居和非兄弟邻居。兄弟邻居是指与当前单元共享同一父单元的相邻单元；非兄弟邻居是指与当前单元不存在兄弟关系的单元。
2.4.1 同级邻居
例如，红色单元 (10,01,01) 的兄弟邻居包括黄色单元 (10,01,00) 和橙色单元 (10,01,11)。这三个单元 Morton 编码的前四位二进制数字 (10,01) 相同，表明它们是兄弟单元。红色单元第五位二进制数字 0 代表 x 方向，第六位 1 代表 y 方向。若要查找 x 方向的邻居，只需保留第五和第六位二进制数字，将 0 取反得到 1，即可找到 x 方向上的兄弟邻居 (10,01,11)，即橙色单元。
若要查找红色单元的非兄弟邻居——绿色单元 (00,11,11)，可在两个方向上进行位运算。例如，确定红色单元 (10,01,01) 沿 x 轴方向非兄弟邻居的过程如下：
确定沿特定轴方向的位翻转层位置。 该单元在最高细化层的 Morton 编码为 01，即 x 方向的二进制位为 0。向上追溯至父层，二进制编码为 01，x 方向二进制位仍为 0。继续向上追溯至祖父层，二进制数为 01，x 方向二进制位变为 1。因此，该层为位翻转层。即从左至右，x 方向的二进制位为 (1,0,0)，第一层为位翻转层。
在对应位翻转层翻转二进制位。 对于红色单元 (10,01,01)，其 x 方向的位翻转层已确定为第一层。为确定红色单元的非兄弟邻居，将红色单元 x 方向的 Morton 编码 (1,0,0) 从第一层起全部取反，得到 (0,1,1)。y 方向的 Morton 编码 (0,1,1) 保持不变，从而得到 (00,11,11)，恰好对应绿色单元。
2.4.2 跨层邻居
若要确定给定方向上细化层数不同的邻居，采用以下方法，以粉色单元 (11,00,00) 在 y 方向上的非兄弟邻居为例：
类似第 2.4.1 节，粉色单元的 Morton 编码为 (11,00,00)，确认位翻转层为第二层。因此，对 y 轴 Morton 编码 (1,0) 取反得到 (0,1)，x 轴保持不变，从而得到 Morton 编码 (10,01)。
位于 (10,01) 的单元细化层数高于当前粉色单元。因此，需在第三层重复最后两位数字，得到 (10,01,01)，对应红色单元。红色单元是粉色单元的跨层邻居之一。若要找到另一个邻居，只需找到红色单元的 x 轴兄弟邻居，即将 01 中的 x 轴位取反，得到橙色单元 (10,01,11)。
由此，粉色单元在 y 轴方向上的两个跨层邻居 (10,01,11) 和 (10,01,01) 均已完整确定。
Morton 编码的优势
数据压缩： Morton 编码的构建方法相对简单，将表示多维数据的坐标值转换为单一整数。通过对 Morton 编码进行位运算，可以便捷地获取网格单元的各类几何和拓扑信息，从而显著减少存储和传输的数据量。
数据局部性： 经 Morton 编码排序后，相邻网格单元往往对应物理上连续的位置。这提高了缓存命中率和内存带宽利用率，从而更好地发挥现代处理器的并行计算能力，执行更复杂的计算任务。
计算效率： Morton 编码排序大幅减少了冗余计算，进而提升了计算速度和整体效率。