命令行运行python_sfepy: python有限元分析模块介绍-【2】从命令行运行命令进行分析...

下面简单介绍一下，采用sfepy进行有限元分析的方法，这里注意几点：

1、sfepy只是进行有限元分析计算，前道工序，几何建模是用其他工具完成，建议用gmesh，后面会撰文介绍，sfepy可以设置材料参数、边界条件等；

2、后道处理是调用mayavi来进行显示，这一块兼容性问题很大，如果不需要，可以不用安装，比如改用matplotlib来显示，或者用pandas，都是很好的替换。

采用 sfepy进行分析的方法，可以参见官网：

Tutorial - SfePy 2019.2+git.02048e0e244c13d4d186df7eabe206906fdfdf84 documentation​sfepy.org

这里简单说明如何从命令行直接运行（把sfepy当作一个黑盒软件 来运行），这里可以用sfepy-run程序运行特定命令，对指定的 problem_description_file进行分析，其语法是：

sfepy-run
usage: sfepy-run [command] [options] [problem_description_file]

Simple wrapper for main SfePy commands.

positional arguments:
{
extractor,
homogen,
phonon,
postproc, （后处理）
probe,
run_tests,
simple  (分析)
}
                      Available SfePy command(s).
options               Additional options passed directly to selected
                      [command].

optional arguments:
-h, --help            show this help message and exit
-v, --version         show program's version number and exit
-w, --window          use alternative (pythonw) interpreter

一般常用的命令就是simple和postproc，分别是分析计算和后处理分析。

这里的problem_description_file是对有限元问题的描述文件，下面以一个示例(即前文中的poisson_short_syntax.py)来说明里面 的基本结构：

这个问题是来描述一个圆柱体(cylinder）在稳定温度场下的热传导分析，适用于Laplace equation来求解，使用到的材料参数只有一个：热扩散系数 。

几何模型文件在data_dir +'/meshes/3d/cylinder.mesh'，描述了一个圆柱体（并划分好了单元） ，其文件结构为：

MeshVersionFormatted 1   (Mesh版本)
Dimension 3              （三维）
Vertices                 （下面是3D点坐标）
354                      （共354个点）
1.0000000000e-01 2.0000000000e-02 -1.2246063538e-18 0
1.0000000000e-01 1.8019377358e-02 8.6776747824e-03 0
1.0000000000e-01 1.2469796037e-02 1.5636629649e-02 0
....
Tetrahedra               (网格划分后的四面体)
1348                     (四面体个数)
29 61 46 30 6            (前面是四面体四个顶点在上面点列表中的序号，最后一个未知)
29 61 58 46 6
29 58 28 46 6
14 71 13 22 6
....

体现的cylinder如下图所示：

下面是poisson_short_syntax.py文件：

from __future__ import absolute_import
from sfepy import data_dir   #这里是获取sfepy安装的目录

#几何模型文件
filename_mesh = data_dir + '/meshes/3d/cylinder.mesh'

#材料参数
materials = {
    'coef' : ({'val' : 1.0},),   #材料的热扩散系数，大概是玻璃
}

#设置积分的区域，以方便后面设置边界条件 
regions = {
    'Omega' : 'all', # or 'cells of group 6'，Omega区域：所有单元
    'Gamma_Left' : ('vertices in (x < 0.00001)', 'facet'),   #左区域Γ，x<0.00001的所有面
    'Gamma_Right' : ('vertices in (x > 0.099999)', 'facet'), #右区域Γ，x>0.099999的所有面
}

#设置变量字段（加在边界条件或初始条件上的负载）的类型
#定义： <name> : (<data_type>, <shape>, <region_name>, <approx_order>)
#<name>: 类型名字
#<data_type>: 变量数值类型，可选： ‘real’ or ‘complex',或 numpy 支持的类型
#<shape>: 负载的维度DOF，可以是：1,2,3
#<region_name>: 这些负载作用在哪个区域，在上面的regions中定义
#<approx_order>: 近似阶数，可以是0, 1, 2, 1B(1 with Bubble)
fields = {
    'temperature' : ('real', 1, 'Omega', 1),
}

#定义变量（用于设置边界条件或负载等）
#定义： <name> : (<kind>, <field_name>, <spec>, [<history>])
#<name>： 名字
#<kind>: 变量种类，可选：'unknown field’(未知字段变量), ‘test field’(因变量) or ‘parameter field’
#<field_name>: 变量字段类型，必须在上面 的fields中定义,
#<spec>: 指定相关，如果变量 是:'unknown filed',需要指定初始值，如果是:'test field'，需要指定相关变量
#<history>: number of time steps to remember prior to current step
variables = {
    't' : ('unknown field', 'temperature', 0),  #变量 t: 自变量，字段类型是上面的'temperature'，也就是
                        #其维度DOF，只有1，也就是要设置或获取其值，要使用't.0'这样的指定
    's' : ('test field',    'temperature', 't'),  #变量s: 因变量,  相当于 s=f(t), 同样只有一维。
}

# Dirichlet 边界条件设置：
#定义：  <name> : (<region_name>, [<times_specification>,]
#              {<dof_specification> : <value>[,
#               <dof_specification> : <value>, ...]
#               } )
#<region_name> : 边界条件加载的区域
#<times_specification>: The times can be given either using a list of tuples (t0, t1) 
#                       making the condition active for t0 <= t < t1, 
#                       or by a name of a function taking the time argument and returning True or False 
#                       depending on whether the condition is active at the given time or not.
#<dof_specification>: 负载的自由度
#<value>: 负载的值
ebcs = {
    't1' : ('Gamma_Left', {'t.0' : 2.0}),
    't2' : ('Gamma_Right', {'t.0' : -2.0}),
}

#定义积分类型和正交规则。
#注意，在这里，积分规则的名字必须以'i'开头
integrals = {
    'i' : 2,  #在3维空间上的二阶正交
}

#设定如何设置组合项：combining term，可以有两个选择：
# Laplace equation, named integral:  dw_laplace.i.Omega( coef.val, s, t ) = 0
# Laplace equation, simplified integral given by order:
#                   dw_laplace.2.Omega( coef.val, s, t ) = 0
equations = {
    'Temperature' : """dw_laplace.i.Omega( coef.val, s, t ) = 0"""
}

#设置线性和非线性的求解器的参数
solvers = {
    'ls' : ('ls.scipy_direct', {}), #线性求解器
    'newton' : ('nls.newton',       #非线性求解器
                {'i_max'      : 1,      
                 'eps_a'      : 1e-10,  #误差范围
    }),
}
#设置求解需要的线性和非线性求解器
options = {
    'nls' : 'newton',
    'ls' : 'ls',
}

上面的文件设置好后，就可以运行sfepy-run进行求解，结果会生成vtk文件，包含所需要的数据。

下一节，我们介绍如何通过ipython的交互式环境，来一步步通过API来设置并进行分析，并通过mayavi来查看不同的postproc的结果。