导读：有限元分析，工业软件里比较重要的基础技术。目前比较流行的comsol，就是基于有限元方法的。那还是3年前的事了，我写了DinoFem，当时是基于1D的，最近有了点时间，补上了2D的有限元程序。期望这次能够一次性做到3D，目前主要还是针对possion方程。

关键字：有限元方法，DinoFem

1、代码结构

2、使用方法及案例

3、结束

代码全部使用了python语言编写的。DinoFem1D和有限元的文章见

【1】跟我学有限元:从inp文件读取网格信息

【2】MESHIO：Python有限元网格文件数据读取与转换。

1、代码结构

代码结构上，我做了一个思维导图，网格剖分，暂时使用了gmsh的网格，和inp文件。实际上，只要是inp文件就可以了。

整个代码结构如下：

目前只是针对线单元来计算的，后期，再慢慢补上其他单元的形式。不知道一晃，是不是又是三年之后了。

年龄越来越大了，就是想沉下心来做一件事。

这三年，虽然说很艰难，实际上不也过来了。

2、使用方法及案例

第一步，生成inp文件

这个目前使用的是gmsh软件，也可以使用其他软件生成inp文件。

第二步，定义系数、外力和边界条件

def solution(x,y):
    return sin(2*pi*x)*cos(2*pi*y)


def ceo_fun(x,y):
    return 1


def force_fun(x, y):
    return 8*pi**2*sin(2*pi*x)*cos(2*pi*y)


def g(x, y):
    return solution(x, y)


def boundary_type_dict():
    a = {}
    a.update({"bottom": -2})
    a.update({"left": -1})
    a.update({"right": -1})
    a.update({"up": -1})
    return a

第三步，求解计算

将上面的inp文件，和函数，并输入一个vtk的输出结构的文件名

def run(inp_file, ceo_fun, force_fun, g,a, result_vtk_filename):
    pt = Mesh2D()
    pt.get_pt_from_inp(inp_file)
    # 设定边界条件
    pt.boundary_type_update(a)
    #  求解器输入网格，方程系数，外力函数和边界函数
    fem = Fem2DLinear(pt=pt,
                      ceo_fun=ceo_fun,  # 方程的系数
                      force_fun=force_fun,  # 外力函数
                      g=g,  # 边界函数
                      )
    fem.solve()
    # 获得解
    v = fem.result
    # 计算误差
    error = [abs(v[i] - solution(pt.P[i][0], pt.P[i][1])) for i in range(pt.number_of_nodes)]
    save_post_vtk(inp_file=inp_file,result=v,error=error,vtk_file=result_vtk_filename)

第四步，可视化

if __name__ == '__main__':
    inp_file = r"D:\Program Files\bConverged\examples\gmshtut\t4.inp"
    result_vtk_filename = "t3_refine4.vtk"
    a=boundary_type_dict()
    # run(inp_file, ceo_fun, force_fun, g, a, result_vtk_filename)
    post_show_from_vtk(result_vtk_filename,"result")

最后就可以去看到结果了。

精确解是：sin(2*pi*x)*cos(2*pi*y)

输出结果：

误差显示：

由于网格量并不多，这个精度还是可以接受的。

只是我们可能要分析一下，为什么在那个区域的误差就那么大。

3、结束

有限元涉及的范围非常广，这几十年来，其实我们本身的有限元人才并不缺，缺的是作为商业软件的商业环境。

期望，我们的有限元软件能做的更好吧。