跟我学有限元：从inp文件读取网格信息

导读：前一篇文章介绍了Meshio的库的使用方法。今天继续讲，如何使用meshio从inp文件读取有限元计算需要的网格信息和边界信息。

inp文件格式，是Abaqus有限元软件使用的格式，开源软件的calculix也是使用了inp文件格式。另外，gmsh网格生成器也提供了生成inp文件的输出工具。这也是我为什么首选了inp格式的原因。

所有的代码，将在github公开，大家可以查找smtmobly/DinoFem

再远的路，都是一步步走出来的。

而，每一步，都很近。

1、有限元网格的基本信息

DinoFem所使用的有限元框架，是基于密苏里科技大学何晓明教授在天元基金的公开课《有限元编程基础》的框架下完成的。为了简单计，我们写程序之初，全部以线性三角单元来进行。

有限元的主要信息包含：

数据信息：    points  存储点的坐标信息    points_size 点的数量信息    cells  存储 单元信息，对应一个点编号列表    cells_size 单元数量信息    boundary_mat 信息，改信息包含3个部分        - 边界单元所在cell的编号        - 边界单元的点的编号的列表    boundary_size 边界上单元的数量函数重定义：    update_boundary_info(self) 设定边界信息矩阵    set_boundary_ibc(self,bc,btype)  设定不同边界上的边界类型的信息

为此，我们写了一个基类DinoMesh

import meshioimport numpy as npclass DinoMesh:    """    这是有限元网格处理的基类。必须包含如下    数据信息：        points  存储点的坐标信息        points_size 点的数量信息        cells  存储 单元信息，对应一个点编号列表        cells_size 单元数量信息        boundary_mat 信息，改信息包含3个部分            - 边界单元所在cell的编号            - 边界单元的点的编号的列表        boundary_size 边界上单元的数量    函数重定义：        update_boundary_info(self) 设定边界信息矩阵        set_boundary_ibc(self,bc,btype)  设定不同边界上的边界类型的信息    1、初始化时只处理网格信息，只设定边界,默认边界值为dirichlet边界条件    2、边界设定在单独的函数set_boundary中进行    3、边界条件设定的值，由类变量__boundary_type字典来设定。        d 代表dirichlet边界        n 代表neumann边界        r 代表robin边界    """    __boundary_type_keys = ['d', 'n', 'r']    __boundary_type_values = [-1, -2, -3]    __boundary_type = dict(zip(__boundary_type_keys, __boundary_type_values))    def __init__(self):        self.data = None        self.faces = []        self.faces_size=len(self.faces)        self.volumes = []        self.volumes_size = len(self.volumes)        self.boundary_nodes = None        self.update_boundary_info()    def update_boundary_info(self):        pass    @property    def boundary_type(self):        return self.__boundary_type    @property    def points(self):        return    @property    def points_size(self):        return    @property    def cells(self):        return self.volumes    @property    def cells_size(self):        return self.volumes_size    @property    def boundary_mat(self):        return self.boundary_nodes    @property    def boundary_size(self):        return self.faces_size    def face_list(self, bc):        return    def set_boundary_ibc(self, bc, btype):        pass

在这个基类中，我们将有限元在网格和边界上的几何信息，都包含在这个基类里了。只是，没有实现。下面我们就以inp格式为案例，进行具体化。

也就是对这个DinoMesh的抽象类进行具体化。

2、Inp2Mesh类

首先初始化的时候，我们需要输入一个文件名，也就是inp所存储的文件。我们的所有的信息都来自于这个inp文件，所以初始化文件的参数，就是这个文件名

    def __init__(self, inp_file):        super().__init__()        self.data = meshio.read(inp_file)

我们需要调用一下父类的初始化函数，父类的初始化函数，其实只是定义了一些变量。

通过meshio库，读取inp文件。

接着，我们继续做一些初始化的工作

        self.faces = self.data.get_cells_type("triangle")        self.faces_size=len(self.faces)        self.volumes = self.data.get_cells_type('tetra')        self.volumes_size = len(self.volumes)        self.boundary_nodes = np.zeros((self.faces_size, 5))        self.update_boundary_info()

（1）把边界上的面的信息提取出来，使用了meshio的get_cells_tyoe()函数。因为是三角元，所以只读取了triangle 类型。

（2）读取体积单元，同样对应的是tetra 四面体

（3）定义边界节点矩阵，这个矩阵包含很多信息，我们后面定义了一个属性boundary_mat直接指向了这个boundary_nodes矩阵

（4）对boundary_nodes矩阵进行填充。

整个步骤就这么多。但是似乎我们都在做边界上的事。

其实，meshio读取信息的时候，points，cells， 这些信息都有了，最麻烦的就是边界信息。

3、填充边界矩阵 （update_boundary_info函数）

    def update_boundary_info(self):        for i in range(self.faces_size):            for j in range(self.volumes_size):                if all([e in self.volumes[j] for e in self.faces[i]]):                    self.boundary_nodes[i, 0] = -1 #类型标识                    self.boundary_nodes[i, 1] = j   # 所属cell编号                    self.boundary_nodes[i, 2:] = self.faces[i]  # 三角元对应的节点的编号列表

这个函数的定义非常简单，也就是对所有的边界上的三角单元寻找他所属于的体积元，也就是cell的编号。同时给他一个类型的标识。

4、设定边界类型 （set_boundary_ibc 函数）

def set_boundary_ibc(self,bc,btype):    value = self.boundary_type[btype]    for i in self.face_list(bc):        self.boundary_nodes[i,0]=value

也就是找到，对应的边界bc，他的面列表对应的boundary_nodes中的边界类型设定值。

bc表示的你为你的边界所起的名字。例如上面的图中，我们给边上一圈命名为face1.

边界类型我们暂时分成了三类：

这三个边界条件是偏微分方程中，常用到的，绝大部分都是由这三种边界条件混合而成。

5、结论

这是一个从inp出发的一个接口程序，将inp网格文件通过meshio和我们的Inp2Mesh类完成了有限元网格信息的输出。

最后贴出全部的，Inp2Mesh的代码

import meshioimport numpy as npfrom DinoFem.fem3d.meshbase import DinoMeshclass Inp2mesh(DinoMesh):    """    1、初始化时只处理网格信息，只设定边界,默认边界值为dirichlet边界条件    2、边界设定在单独的函数set_boundary中进行    3、边界条件设定的值，由类变量__boundary_type字典来设定。        d 代表dirichlet边界        n 代表neumann边界        r 代表robin边界    """    def __init__(self, inp_file):        super().__init__()        self.data = meshio.read(inp_file)        self.faces = self.data.get_cells_type("triangle")        self.faces_size=len(self.faces)        self.volumes = self.data.get_cells_type('tetra')        self.volumes_size = len(self.volumes)        self.boundary_nodes = np.zeros((self.faces_size, 5))        self.update_boundary_info()    def update_boundary_info(self):        for i in range(self.faces_size):            for j in range(self.volumes_size):                if all([e in self.volumes[j] for e in self.faces[i]]):                    self.boundary_nodes[i, 0] = -1                    self.boundary_nodes[i, 1] = j                    self.boundary_nodes[i, 2:] = self.faces[i]    @property    def points(self):        return self.data.points    @property    def points_size(self):        return len(self.data.points)    @property    def cells(self):        return self.volumes    @property    def cells_size(self):        return self.volumes_size    @property    def boundary(self):        return  self.boundary_nodes    @property    def boundary_element_size(self):        return self.faces_size    def face_list(self, bc):        return self.data.cell_sets_dict[bc]['triangle']    def set_boundary_ibc(self,bc,btype):        value = self.boundary_type[btype]        for i in self.face_list(bc):            self.boundary_nodes[i,0]=valueif __name__ == '__main__':    a=Inp2mesh("../resources/fem_test003.inp")    a.set_boundary_ibc('face1','n')    a.set_boundary_ibc('face2', 'r')    print(a.boundary_nodes)    print(a.boundary_nodes[a.boundary_nodes[:, 0] == -1, :])

来源：张麟博士