非结构网格下寻找n层边缘网格单元的算法

目的：在一张非结构网格当中寻找最外n层的边缘网格

原理：在二维情况下，如果有两个网格相连，那么一定会在中间有一条公共边，这条公共边属于这两个网格，但是对于边缘网格，一定存在一条边只属于该网格自己。三维情况下，公共边对应为公共面，即两个相邻网格单元中间存在一公共面，该公共面同时属于这两个网格单元，而对于边缘网格单元，必定存在一个单独的面只属于一个网格单元。通过遍历每个网格，统计每一条边（面）的遍历次数，筛选出来只遍历一次的边（面）就可以找到边缘网格。

概述：非结构网格与结构化网格相对，结构化网格是指每一个网格周围相连的网格网格数量相同，相对的，非结构网格是指每一个网格单元附近有任意数量的网格相连接。非结构网格具有剖分灵活，易于刻画复杂界面和复杂几何的特性。在计算流体力学，固体力学当中都有相当广泛的应用。边界层由于其特殊的性质，一直以来都是研究的热点，该算法可以将一套非结构/结构网格当中的边缘n层网格提取出来。

（图片摘自百度图片）

实现语言：Python

依赖库：numpy

数据结构：每个网格单元存储（面/边）对应的编号，例如两个四边形的四条边的序号分别为0，1，2，3 和 3，4，5，6。那么单元数组的格式如下：

element = [[0,1,2,3],[3,4,5,6]]

1. 从寻找的一层边缘网格说起

    正如前面原理所述，所有最外层网格必存在至少一条边（面）是非公共边（面），也就是说该边（面）的序号只会在整个element数组里面出现一次，其余的公共边（面）会出现两次。那么只需要寻找在该数组中只出现1次的数并保存下来即可以得到这些单独的边（面）的序号，得到这些边的序号以后再反查他们属于哪一个单元即可以得到最外层网格的单元序号。

    边（面） = {非公共边（面），在element数组中只出现一次} ∪ {公共边（面），在element数组中出现两次}

边缘网格 = {存在至少一条边为（面）为非公共边（面）的网格}  

    假如存在两个四边形网格，element=[[0, 1, 2, 3], [3, 4, 5, 6]]，可以看到3为公共边，它在element里面出现两次，而边 [0, 1, 2, 4, 5, 6]都是非公共边，它们在整个element里面出现两次。它们来自于网格0和1，就是说这两个网格单元都是边缘网格。

程序实现：

    1.1 寻找非公共边（面同理）

        先将element数组展平，然后得到所有的边的编号

element_flatten = element.flatten()

        然后使用np.unique函数去掉重复值，得到所有公共边和非公共边的索引，这里需要使用return_index=True，来返回索引。return_index=True会返回某个数第一次出现的位置，之后若还出现该数字则不予记录。所以该索引记录了非公共边的索引和公共边一半的索引。

edgeList, indexes = np.unique(element_flatten, return_index=True)

        求取上述该索引的补集，就可以得到所有公共边的编号，使用np.delete从element_flatten去除edgeList即可得到补集了。

couple_edges = np.delete(element_flatten, indexes)

        将所有边的序号(edgeList)和所有公共边的序号(couple_edges)做差集，就可以得到所有非公共边的集合。set类函数提供了difference来求差集，set(a).difference(set(b))返回a当中存在而b当中不存在的集合。

lonely_edges = np.array(list(set(edgeList).difference(set(couple_edges))))

        得到所有非公共边之后，遍历这些非公共边出现在哪些网格单元内得到边缘网格的索引。下面代码中pml_indices就是边缘网格的索引。

        pml_indices = []        #pml_indices为边缘网格的索引
        for i, edges_per_element in list(enumerate(element)):
            for edge in edges_per_element:
                if edge in lonely_edges:
                    pml_indices.append(i)

        至此就得到了所有边缘网格的索引。

2. 向n层边缘网格进发

    寻找n层网格和寻找1层网格的原理完全相同，只不过需要在寻找完1层边缘网格以后，将寻找得到的边缘网格从所有网格当中删除，就像剥洋葱一样。这样再在留下来的网格当中去寻找最外层的网格，如此循环n次就可以最终得到n层边缘网格的目标。在具体实现的时候不能直接使用上面代码的方式得到索引，因为剥去外层网格以后，pml_indices对应的索引是改变之后element的索引。这时候需要提前将element备份，通过直接比对单元的边编号来确定它在全局element的位置。

        具体请见下面代码：

def findPML(n_layer=1, element=None):
    import numpy as np
    element_bak = element
    element = np.array(element)
    pml_indices = []        #边缘网格的索引
    pml_elements = []        #边缘网格的网格单元
    for i_layer in range(1, n_layer+1):
        element_flatten = element.flatten()
        pml_indices_once = np.array([])        #第i_layer层在剥离(i_layer)层边缘网格后的网格索引
        edgeList, indexes = np.unique(element_flatten, return_index=True)
        couple_edges = np.delete(element_flatten, indexes)
        lonely_edges = np.array(list(set(edgeList).difference(set(couple_edges))))
        for i, edges_per_element in list(enumerate(element)):
            for edge in edges_per_element:
                if edge in lonely_edges:
                    pml_indices_once = np.append(pml_indices_once, i)
        pml_indices_once = np.array(np.unique(pml_indices_once), np.uint8)
        pml_elements.append(element[pml_indices_once].tolist())
        element = np.delete(element, pml_indices_once, axis=0)
        pml_indices_once_list=[]
        for _ in pml_elements[i_layer-1]:
            pml_indices_once_list.append(element_bak.index(_))
        pml_indices.append(pml_indices_once_list)
        if len(element) == 0:        #如果剥离i_layer层以后不存在网格了，那么break提前退出循环
            print("delete %d layer(s)" % i_layer)
            break
    return pml_indices, np.array(pml_elements)

if __name__ == '__main__':
    element = [[0,3,4,7],[1,4,5,8],[2,5,6,9],[7,10,11,14],[8,11,12,15],[9,12,13,16],[14,17,18,21],[15,18,19,22],[16,19,20,23]]
    indices, elements = findPML(n_layer=3, element=element)
    print(elements)
    print(indices)

       

        得到的结果如下：

[[[ 0  3  4  7]
  [ 1  4  5  8]
  [ 2  5  6  9]
  [ 7 10 11 14]
  [ 9 12 13 16]
  [14 17 18 21]
  [15 18 19 22]
  [16 19 20 23]]]
[[0, 1, 2, 3, 5, 6, 7, 8]]

        将n_layer改为3:

delete 2 layer(s)
[ [[0, 3, 4, 7], [1, 4, 5, 8], [2, 5, 6, 9], [7, 10, 11, 14], [9, 12, 13, 16], [14, 17, 18, 21], [15, 18, 19, 22], [16, 19, 20, 23]]
 [[8, 11, 12, 15]]]
[[0, 1, 2, 3, 5, 6, 7, 8], [4]]

        可以看到只剥离了两层以后就退出循环了，返回的数组正确的返回了每一层网格的信息。