怎样开发一款有限元软件，从哪些方面学习？

开发有限元软件是一个复杂而又艰巨的工程，今天我们来看看知乎上的大牛都是怎么看的。

行业大牛ZifengYuan认为：

这是一个很大的问题。我自己从零写过一个功能还算详全都有限元代码，借着自己的一些经验，尽可能回答一二。

动机

市面上，已有的著名商业有限元软件占据了大部分工业界的份额；学术界可以选择商用软件（当然可以选择教育版许可证），亦可以选择一些常见的开源有限元代码做二次开发。基于以上，大多数问题都可以用金钱和/或时间解决。但为什么还有许多人源源不断的选择自己开发有限元软件？

从零开始设计有限元软件的动机或者说好处简单来说有如下几类：

1）版权问题。正版软件哪怕是教育版的许可证也是很贵的，而且如果教育版的功能还有限制，是很制约研究水平的。譬如ABAQUS教育版不能使用UMAT等API功能，对自由度也有上限规定。另一方面，若自己开发商用上层代码，不能随意使用别人的开源底层有限元程序，因为这也可能是有开源协议所保护的。

2）研究方向的问题。做计算力学研究的话，我们经常需要开发诸如新的单元技术，本构模型等。一些新的计算技术可能无法在商业软件二次开发中实现。而大型的商业软件可以店大欺客，也就是他们并不会因为你们的学术需求定制额外的功能。另外有一些研究方向，譬如反问题、随机问题、优化问题等，需要大规模反复对同一个模型进行有限元计算，这时候若是底层有限元代码是自己写的，无论是建模还是后处理将大大减少工作量。

3）想通过编程更好理解有限元思想。有限元作为一种数值求解偏微分方程的计算方法，书本的内容再详细，最后实现还是通过程序。通过编写有限元程序，可以更好的理解其思想，也能在解决实际问题中更好的发现问题解决问题。

4）好的有限元程序可以加速整个团队的研究。我们组从2015年开始基于我的代码做了一些工作，一共发表了5篇论文（两篇IJNME，一篇FINEL，一篇JAM，一篇IJF），以及两次会议报告。后续还有一些工作正在研究中。我在其中的作用基本是顾问，给予一些简单的代码语法框架等指导，然后同组人员做自己的专门的研究。

语言

使用什么样的编程语言编写有限元程序是一个很有意义的话题。一开始是毫无悬念的FORTRAN 77（FEAP，Warp3D）；后来随着面向对象概念的兴起，出现了不少C/C++编写的代码（OOF3D，dealii，OOFEM）。FORTRAN 77格式残留了不少，因此后来的开发者被迫向其兼容而痛苦不堪。

从优化的角度说，程序运算的速度其实和机器有关，而原则上大多数语言都可以优化到同一个水平。因此，考虑语言的时候，计算效率原则上不应该成为取舍的主要原因。需要考虑的因素主要在于对语言的熟练程度，便于维护，便于合作，能跳过语言的一些坑，能熟练debug。

我本人一直使用的是Intel Fortran，现在已经升级到15.X，已经涵盖大部分Fortran 2003/2008标准。如今的Fortran和当初的FORTRAN 77已经大不相同，也有面向对象的功能。当然一个语言不同的人也会演变成不同的风格。我个人很不喜欢FORTRAN 77， 但也反对对Fortran语言的全盘否定，因为自打Fortran 90开始，这门语言正在努力的更新换代。尤其从Fortran 2003/2008开始，Fortran以及支持面向对象编程，尽管在有限元代码中，面向对象的不少特性我们用的不多。

我自己从初一开始参加OI，从Pascal开始学习，拿过一些省一等；后来在本科时候学了点C，C++完全是自学，基本不会面向对象；出国后抱着大腿参加过一次ACM-ICPC World Final（哈尔滨那次）。十几年的非CS专业的编程历史是一段很有意思的经历。知乎上CS方面的话题基本上看不懂；刚开始编写有限元软件的时候，不懂维护，不懂版本控制，不懂架构。甚至我编写的环境是Win+VS——大写的”土法炼钢“不是么。后来开始慢慢充电，看软件工程方面的书，思考架构，思考怎么维护代码。于是软件的可读性、可维护性、可拓展性才慢慢提高。

开始

万事开头难。在刚开始写有限元代码的时候，面对空白的project，该如何开始？不妨先思考一下这些问题：

0）编程语言。

1）是针对什么问题的有限元？

有限元方法适用于不少偏微分方程的数值解法。结构力学可以用，连续介质力学也可以，而对于物质的扩散和热传导问题，有限元更是具有有限差分无法比拟的优点。

2）是一维、二维，还是三维问题？

关系到前后处理、单元库的编写。

3）对于有限元问题，如何定义合适的类（class）？

有限元里面的概念纷繁复杂。数学上的公式更多的是算法范畴，而数据结构那块需要合适的将信息分割成合适的类。

4）如何规定输入文件的格式，如何编写一个简单的parser？

我们编写通用有限元代码不是为了某一个特定模型，而是某种意义上通用问题。因此，我们需要定义合适的输入文件的格式：简单，但又有一定的语法规则。

5）用什么前后处理软件？

基本上我们不会去编写前后处理的部分，因此我们需要选择一款合适的软件来做网格划分以及后处理显示。

类的设计-0

这一篇主要讲有限元程序中关于数据结构的设计。环境是集成Fortran 2003/2008（90/95也适用）的Intel Fortran。

在Fortran 2003/2008中，类（class）及其类成员的继承、覆盖等简单的面向对象概念终于被支持。在我自己设计的程序中，面向对象的概念并没有广泛使用，最主要的是我并没有接受过很好的面向对象程序设计的训练，熟练度差；其次是我们经过简单测试，Fortran的面向对象技术会降低计算效率。因此，尽管我这里称之为类，但实际上仍然是90/95环境下的derived data type。

和C++相比，Fortran缺乏头文件（header）的引用模式，而是采用module的引用（注意，还不能相互引用）来模块化代码（类似C++的namespace）。因此，Fortran里面通过derived data type来设计一个类其实是很不爽快的。举个例子，我们有两个类A和B，A包含一些操作需要B类的信息，而B包含一些操作需要A类的信息。在C++里面就很好实现，此处不多叙述。而在Fortran里面就不能很直接的实现。譬如这样做就是非法的：

MODULE A_subs_mod

  USE B_subs_mod

IMPLICIT NONE

  TYPE :: A_class

    ! members of A class

  END TYPE A_class

CONTAINS

  ...

END MODULE A_subs_mod

MODULE B_subs_mod

  USE A_subs_mod

IMPLICIT NONE

  TYPE :: B_class

    ! members of B class

  END TYPE B_class

CONTAINS

  ...

END MODULE B_subs_mod

[注：可能2008标准解决了一些，但当下不少编译器还没有完全涵盖2008标准]

因此，折中一下，数据结构和程序文件如下设计，对于任意一个类A，我们设计A_def.f90, A_cmn.f90, A_mod.f90等一系列程序文件，其中A_def.f90包含最简单的类的定义模块：

MODULE A_def

IMPLICIT NONE

  TYPE :: A_class

    ! members of A class

  END TYPE A_class

END MODULE A_def

由于Fortran不支持static成员，也没有枚举类型，因此，额外的常量参数也可以写在这个文件。

A_cmn.f90包含一些简单的关于类A的数据操作，譬如成员的赋值，指针的allocation, dellocation, association, 和nullify，可变数组（allocatable array）的allocation, 和deallocation等：

MODULE A_cmn

  USE A_def

IMPLICIT NONE

CONTAINS

  SUBROUTINE A_set_A_member1(obj, A_member1_value)

      TYPE(A_class),INTENT(INOUT)    :: obj

      DATATYPE(KIND=?),INTENT(IN)    :: A_member1_value

      A%member1 = A_member1_value

      RETURN

  END SUBROUTINE A_set_A_member1

END MODULE A_cmn

最后的A_mod.f90包含一些高级一些的subroutine，可能需要其他的类B作为input argument。但注意这里我们只允许引用B_def或者B_cmn，而一般避免引用B_mod，从而规避交叉引用。若是类A和类B是包含与被包含关系，则允许引用B_mod。

这里的基本程序文件譬如A_def.f90, A_cmn.f90等程序可以通过程序自动生成。我自己的写的程序生成器可以生成大约10类程序文件供使用，以及可以生成合适的帮助文档供查阅。程序自动生成可以避免bug的产生，也大量节省了编程时间。

类的设计-结构

有限元程序中类的结构大致分为如下几块：

1）IO：包含文件（file），语法分析（parser）等；

2）数学：包含张量（tensor），amplitude，orientation等；

3）库（library）：包含材料库，单元库，物理场库。库的概念以后会详细讲。

4）有限元数据：包含网格（mesh），材料，物理场，分析步，后处理请求等。

以后几个篇幅简单介绍一下3）和4）的设计。

类的设计1-关于网格（mesh）的类的设计

网格是有限元程序中非常重要的信息来源。网格不仅仅定义了几何信息，譬如结点的坐标以及单元结点编号，在定义材料（及材料坐标架）、边界条件、载荷以及后处理的时候，都需要和网格建立一定的关联。因此，在程序中，如何设计网格的类及其子类，是具有很大意义的。

下面的讨论基于多物理场和几何非线性的有限元代码设计。

关于网格的基本知识

我们不妨从介绍网格的基本知识入手。通常来说，网格包含结点（node）和单元（element）。结点有编号以及坐标信息；而单元有单元编号、单元结点编号和单元类型等信息。拥有相同结点数的单元可能类型是不同的，譬如常见的完全积分的8结点六面体单元C3D8与缩减积分的8结点六面体单元C3D8R。

深入一点说，结点在不同物理场拥有不同意义的自由度（或者说场变量field variable），譬如机械场下，每个结点有三个位移自由度；在温度场下，有一个温度自由度；在扩散场下，有一个物质的量自由度。而对于给定的单元类型，我们也就知道了其积分点的性质——譬如积分点的数量、在参数坐标下的位置、权重等等。

在有限元中，整体刚度矩阵对于隐式求解来说非常重要，整体刚度矩阵的尺寸不考虑边界条件的情况下等于所有自由度之和，一般来说是非常巨大的。这也就是为什么我们常常用稀疏矩阵来存储整体刚度矩阵，也就是存储非零元素的位置和值。因此，我们需要在网格这个层面寻找一些信息，用于定义稀疏矩阵的sparsity关系——事实上，对于任意的结点，我们需要找到所谓的其“相邻”结点。两个结点是相邻的，如果这两个结点同时出现在一个单元内。在整体刚度矩阵中，相邻结点自由度所决定的位置就是非零元素。

网格另一个很重要的信息就是集 合。集 合可以是结点的集 合，也可以是单元的集 合。对于线载荷和面载荷来说，我们还分别需要线单元和面单元的集 合。这是因为我们得通过积分来计算线载荷和面载荷所对应的等效结点力，因此，光给出载荷所在的结点集 合是不够的。

类的设计1-1：结点

结点相对来说最简单。编号、坐标、相邻节点编号数组以及包含这个结点的单元编号数组。

编号和坐标不说，基本上前处理软件都能给出。相邻节点编号数组用于计算整体刚度矩阵非零元素的位置。而最后一个，包含这个结点的单元编号数组可以用于后处理的应力磨平计算中。

类的设计1-2：单元

单元类稍许复杂一点，包含编号、单元结点编号数组作为基本信息，以及单元体积和特征长度等几何信息。此外，为了支持laminate单元，我们额外提供ply的厚度、方向，以及stacking direction信息。

类的设计1-3：积分点

在程序中，积分点的信息设计成与单元分离。这是为了应付可能出现的复杂情况。譬如continuum单元是Gaussian积分点，而cohesive单元是Newton-Cotes积分点，含义不同，信息也不同。与单元分离存储利于更好的管理数据。

类的设计1-4：集 合

这里的集 合很大程度上依赖于Gmsh里面的信息。在Gmsh这个前后处理软件中，几何模型分为点、线、面和体等entity，在网格定义中，Gmsh非常贴心的给出所有entity的mesh信息。譬如说，一个box有8个点，12条边，6个面，1个体，那么在Gmsh的网格文件中，就有8个结点，以及8个含有一个点的单元，12组一维单元，6组二维单元，和一组三维单元。因此，集 合这个类，就是Gmsh中一个entity对应的信息。而随后的材料与单元的关系、边界条件和载荷的定义、后处理的定义，都只与entity的维度和编号相关联，而无需定义繁琐的结点或者单元列表。

科学匠人：

关于语言的问题，补充一下先前我不太推荐c++主要因为复杂度高，对于开发人员要求较高。同时要求开发人员既是计算力学专家，又是c++高手，这对hr来说是艰巨的任务。不过随着技术发展，c++在一些领域有很强的优势体现出来，那是非用不可了。

------------本人写过cfd程序，虽然和fem不是一码事，但是从计算力学求解器的角度看本质也差不多。两条路：1 从头自己写。2 在开源库基础上开发。前者适合快速实现功能，但是如果经验不足会花费很多时间修改重构。后者适合你有足够的时间尝试和学习，因为那么多开源库并非都名符其实，或者适合你的需要。

对于从头开发，我的经验是：

1 循序渐进，先考虑串行单块。设计好松耦合的模块：几何处理，边界条件，单元内rhs和lhs构造，全局lhs构造，线性方程组求解.... 同时设计好层次结构，形成高层模块基于低层模块的结构。单块串行的设计非常基础也非常重要。在这个阶段设计决定了整体设计的好坏。而且一个简单的单块版本适合于一些新方法和新思路的快速尝试。

2 构造较为独立的并行框架，结合第一步的单块求解器，形成另一个独立的并行求解器。

3 基于各个模块，可以采用python脚本编写高层次逻辑，实现灵活的定制应用开发。这似乎可以成为一种DSL了，不过这并不是必须的。

4 充分利用现有的基础库，线性代数方面有blas/lapack，petsc，superlu等，输入输出方面有xml，json以及hdf5或者cgns，分区方面有metis，还有日志、trace等等

5 版本控制是必须的，git目前是标准了，svn也不错

6 前处理我没有搞过，都是用现有的工具生成的网格；后处理也是如此。

7 准备几个常用测试算例，每次修改后都必须通过全部测试。

8 语言的选择。这个最容易引起争议，因为这是最容易喷的话题。我的建议是fortran95或者c这样成熟且简单的语言，而不是fortran2003和c++这种新的还未成熟稳定或者复杂容易出错的语言。

------------------

补充：文档工作也非常重要，推荐doxy工具，直接从代码生成文档。

蚍蜉撼大树：

开发一款有限元软件基于两个方面。

首先，需要具备专业的有限元知识与理论，就是我们所说的算法思想。比如理解位移法，薄板理论，梁理论，等等。当然，也可以在想编写的软件内加入最新研究成果，比如界面粘结等先进理论等等，成为国内首创应用软件，就更具有诱惑力~~

然后，你需要具备一门或一门以上的程序设计语言。传统的有限元程序采用的是面向过程的Fortran程序，当然，近年来面向对象的C++语言更多的应用于有限元软件设计内。当然，如果你对其他程序语言感兴趣，比如C语言，Java都可以编写出优秀的应用软件。甚至你完全可以基于JavaEE用Android开发出相应的有限元软件。

软件的开发最主要的在于算法的思想，其次，在于界面的显示。对于日臻成熟的算法思想下，现今都在比拼界面的优化。提倡界面的简洁，易用性。

题主还说对于目前的某些有限元软件如何做进一步的改进，一般而言是以插件的形式改进。

yanfeng39：

这是一项非常复杂的工程。。。如果你想做成像国外老牌厂家那样的软件，那么凭一个人的能力基本不可能；即使是一些优秀国产软件，背后也是一个公司在支撑。如果真要个人开发软件的话，那么再说一说：软件应该具备前处理、核心算法模块及后处理，并且易于扩展。简单的前处理包括读入数据文件，复杂一点的前处理则包含显示、交互建模、网格划分功能，这些可以用OpenGL及开源程序库实现，但要做到无所不包则不太可能。核心算法则与有限元相关，要尽量做到通用、可扩展，例如扩展材料库、单元库等等。后处理包括结果文件输出、位移云图／应力云图显示、时程曲线显示等等。

可以看出，要做好一个有限元软件是十分费时费力的。