跟我学有限元:从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 meshio
import numpy as np


class 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.

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

• Dirichlet边界 对应的类型值 -1
• Neumman边界 对应的类型值 -2
• Robin 边界 对应的类型 -3

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

5、结论

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

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

import meshio
import numpy as np
from DinoFem.fem3d.meshbase import DinoMesh

class 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]=value


if __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, :])