by HPCer @ Beijing，2020-6-21

说在前面的话

光说不练假把式，这次来点干货，我们实战一次。

7年前还在读博士的时候，在科学网写过一篇博客，用Matlab以不超过150行代码实现了一个简单的二维弹性有限元程序，我们不妨拿这个小例子开刀（Matlab代码可以从 http://blog.sciencenet.cn/blog-3760-678647.html 下载）。

我这个小程序没法在Octave下直接运行，因为涉及到了PDE Toolbox的使用。很多人说Matlab不可替代的是Toolbox工具箱，我也想在这篇文章argue一下，不是所有的工具箱都不能替代的，其实作为Matlab的替代品，Octave也仿照Matlab构建了大量的package，大部分是第三方作者开源贡献的，不信您上 https://octave.sourceforge.io/packages.php 瞧瞧，是不是有很多相似的面孔。我想借用毛主席的名言，“没有调查就没有发言权” （：

Matlab PDE工具箱的替代

我们今天为了完成这个小迁移就来找找PDE Toolbox的替代品。

Octave用来替代PDE Toolbox的工具包主要有两个，一个是生成网格用的msh（https://octave.sourceforge.io/msh/index.html），一个是显示网格和计算结果用的fpl（https://octave.sourceforge.io/fpl/index.html）。

先来看看生成网格的msh。

msh生成非结构化网格实际是后台调用的gmsh，刚好我的macbook上已经有了，如果没有该程序可以简单的用homebrew安装一下，执行：

brew install gmsh

安装工具包的过程很简单，将octave目录转移到下载的工具包所在目录后，执行pkg install xxx，我们这儿执行：

pkg install splines-1.3.3.tar.gz

另外msh工具包还要依赖于spline样条函数工具包，也需要提前下载好，同样执行：

pkg install msh-1.0.10.tar

fpl工具包安装也类似，同样执行：

pkg install fpl-1.3.5.tar.gz

至此，完成我们所有准备工作。

Matlab程序迁移

一个不幸的消息是Matlab工具包PDE Toolbox中的initmesh函数没有直接的替代品，需要使用msh工具包中的函数改写。

fpl的官方文档即给出了和我们Matlab程序生成三角形网格类似的例子，直接抄过来，如下图所示：

Octave改写的网格生成程序

熟悉gmsh的同学可以发现，这个程序就是简单的调用了一下gmsh而已，连命令都一样(: 

执行该程序可以得到我们之前在Matlab中一样的网格，数据结构也是和Matlab PDE工具箱中的一样，结果如下图所示：

Octave生成的二维网格

好了，Matlab中的initmesh函数用我们上面写的程序换换就好。plf工具包生成的网格都在mesh结构体变量里，我们把它简单的赋值出来给p、t、e就好，这也是我们这个小程序移植过程中唯一的改动。

除此以外，还有个需要注意的地方，Octave里第三方工具箱使用前需要pkg load一下。

好了，我们把所有的改动都写在一起：

Matlab程序的改写

运行一下试试，然后我们画图：

Octave画图代码 Octave画图结果

这就完成了Matlab程序的迁移了，是不是很简单。

所以说，Matlab被禁大家不要怕，迁移很容易，如果您不会，那就找我呀（：

用FELAC写一个二维弹性有限元程序

最后我还想说，其实写一个弹性力学有限元程序还有更简单的方法 - 

这就是用我们有限元语言及其编译器FELAC软件。

管你MATLAB是不是被禁，管你OCTAVE是不是有相关的工具箱，用FELAC写有限元程序真的只要不到一分钟 ^-^

FELAC公式库里的FELAC二维线弹性有限元语言脚本 FELAC生成的二维线弹性有限元C语言主程序 FELAC计算结果，x方向位移 FELAC计算结果，y方向位移

附录：增加了一些其他商业软件的计算
ANSYS结果：

ANSYS MECHANICAL计算结果：ux
ANSYS MECHANICAL计算结果：uy

COMSOL结果

COMSOL计算结果：u
COMSOL计算结果：v