获取非结构化多面体网格的相邻单元格图

计算科学 有限体积 开放式泡沫 非结构化网格 网格遍历
2021-12-03 19:27:01

我正在做一个关于在固体立方体上使用有限体积法求解热方程的小项目,我将立方体的多面体网格转换为 OpenFOAM 网格。

我有一个 Python 代码,我在其中解析 OF 案例网格的点、面、所有者、邻居和边界文件,我设法创建了一个哈希列表(让我们称之为CellFacesMap)以将单元格编号映射到其拥有的列表面对并面对它的邻居。现在,我正在尝试生成一个哈希列表(我们称它为CellAdjacentsMap将每个单元与其所有相邻单元映射在一起,以便将离散热方程应用于每个单元(以及生成体积加权因子)。

我现在能想到的唯一方法是访问其中的每个单元格,CellFacesMap对于拥有的面孔,我得到与该面孔相邻的单元格(相邻单元格),对于与同一单元格相邻的面孔,我得到拥有这些面孔的单元格(也相邻的单元格),这意味着我CellFacesMap一遍又一遍地访问每个单元格以检查其面孔。

但是这种方法成本很高(时间方面,为 4000+ 个网格生成地图需要 3.5 分钟)。那么有没有更快的生成方式CellAdjacentsMap呢?

2个回答

使用哈希映射会为所有访问增加 log(n) 复杂度(然后遍历整个网格通常会花费 n log(n)),因此显然它不是最佳解决方案。

现在您的问题是如何有效地存储和遍历邻接信息,即每个单元格的相邻单元格列表。您可以将一个单元格映射到其邻居列表的第一件事:您可以使用数组而不是哈希,这将消除 log(n) 成本。

然后,使用数组,这意味着您需要存储不同长度的列表数组(如果您的网格有任意多面体),这通常不是很有效:每个列表都有自己的内存开销要添加到使用的内存中通过它包含的元素。存储此类信息(即不同长度的列表数组)的一种更有效的方法是“压缩行存储”(称为 CRS、CSR 或耶鲁格式)表示 [1],用于存储稀疏矩阵。

编辑如果您的输入中不存在邻接信息,那么您将需要重新创建它。而不是使用其他答案中建议的散列或关联映射,另一种方法是将所有方面存储在一个数组中然后排序这个数组。排序后,与相邻单元格对应的分面对在数组中是连续的。我观察到它的效率要高得多(至少在 C++ 中,我没有在 Python 中测试过)。另请参阅我对这个问题的回答 [2]。为了使我的答案适应您的多面体网格,构面记录将具有一个构面索引,并且作为键,三个顶点索引(构面中的最低顶点 id 和通过边缘连接到它的两个顶点,按递增顺序)。通过按字典顺序对构面记录进行排序,您将检索数组中连续的一对相邻构面。

算法:

For each cell c
   For each facet f of cell c
       Find v1 the smallest vertex index in facet f
       Find v2 and v3, the vertices before and after v1 in facet f
       If v3 < v2, swap v2 and v3
       Append a new record [v1, v2, v3, c] to the array

Sort the record array with the lexicographic order on (v1,v2,v3)

In the record array, find (by scanning the array) the contiguous 
pairs of record (v1,v2,v3,c1);(v1,v2,v3,c2) with the same (v1,v2,v3). 
For each such pair, c1 and c2 are adjacent.

如果输入数据具有 facet id,则记录只是 (f,c) 对,它们按 f 排序(而不是 (v1,v2,v3))。

[1] https://en.wikipedia.org/wiki/Sparse_matrix

[2]生成元素边数的网格选项(tetgen-triangle)

我假设您从单元定义列表开始,例如,作为定义单元的顶点列表和定义每个单元的拓扑的“类型”。作为每个单元的拓扑定义的一部分,您可以轻松地获得该单元定义的面,例如,作为顶点列表。

创建您需要的 CellAdjacentsMap 的第一步是首先创建一个我称之为 FaceCellMap 的地图,其中列出了每个面所属的单元格(一个或两个)。为此,您可以遍历所有单元格,将每个单元格的面逐个添加到地图中。

形成这个 FaceCellMap 有一些实现细节。首先是你使用什么数据结构?如果您具有通用的 HashMap 功能,则键可以是面定义,值可以是长度为 2 的单元格 ID 数组(其中一个条目可能为零)。第二个问题是你如何比较面孔?您可以使用顶点列表,但您已经实现了比较函数,它对“起始”顶点以及是否以顺时针或逆时针顺序遍历顶点列表不敏感。

但是在创建 FaceCellMap 之后,创建 CellAdjacentsMap 就很简单了。您遍历所有单元格。对于每个单元格,您将遍历其所有面。对于每个面,检查该面的 FaceCellMap 以查看它是否具有相邻的单元格。如果是这样,请将其添加到当前单元格的邻接列表中。