在上篇文章 OpenSees 不同解释器的性能比较中，比较过不同 OpenSees 调用方法的速度不同。直接使用 C++ 调用的速度最快， TCL 调用次之，而用 Python 调用的速度大大减慢。然而， Python 作为当前流行的“胶水”语言，应用非常广泛，并且有很强的扩展功能。因此，使用 Python 快速地调用 OpenSees 是很必要的。

本文提出的加速方法，是先使用 C++ 建立分析模型，并封装为一个函数。然后通过 Python 来调用这一函数，完成分析。

问题定义

本文以一个参数化分析的问题为例。该问题与前文相同。是最常见的 OpenSees 例子，桁架的静力分析。如图所示。

OpenSees性能比较

这一问题中采用的是英制单位，不过不必担心，我们这里不讨论这些数值有多大，只是探索其性能，所以仍保留使用英制单位。

参数化分析中，把水平力从 1 变化到 1,000，再把 A1 的面积从 1 变化到 10，取变化步长为 1 ，则构建了 10,000 个不同的模型。然后我们把这 10,000 个不同模型的顶点 4 的位移进行输出。

使用环境

操作系统

本文使用的环境为 Linux 环境。作者的 Linux 环境为 Alpine 3.12，在 docker 中运行。用户可以自己选择合适的运行环境，MacOS，Windows 都同样支持，只是做法有一定差异。对于差异之处会在文中提到，但是实际情况可能会有所不同。

Windows 用户可以使用 Mingw-w64 来编译 C++ ，也可以选择使用 WSL2 (Windows Subsystem for Linux 2) 来使用 Linux 系统。 Mac 用户和 Linux 用户的编译方法基本相同。

使用 docker

使用 docker 只是为了建立 Linux 虚拟机，并不是必要的。如已有 Linux 环境，可以跳过本节。

对于各平台使用者，都可以使用 docker 来创建虚拟 Linux 环境。具体操作方法，可以参考 docker 的官方文档。本文也作一个简要的介绍。零基础者也可以建立自己的 docker 系统。

首先，按官方文档的方法，安装 docker 。在本地启动 docker server 。打开终端，输入

1
docker --version

如果打印出了版本信息，则安装成功。

下面我们建立一个虚拟 Alpine 3.12 容器，并运行。

1
docker run -it alpine:3.12

run 表示运行一个容器，容器就是一个虚拟的 Linux 实例。

-it 表示交互式。如果没有这两个 flag ，容器会自动退出。加上它，会进入容器中，与容器交互。

alpine:3.12 是镜像。3.12 是版本号。

由于当前的系统中还没有 alpine:3.12 这个镜像，docker 会自动从源拉取，然后再运行容器。

Alpine 是一个只有 5.2M 的 Linux 操作系统，只包含最基础的架构，很多包要自己安装。但是安装方法也很简单。

安装软件包

下面安装所需软件包。这里以 Alpine 的 apk 为例。其它 Linux 系统可以使有对应的包管理器。

在容器中输入

1
2
apk update && apk upgrade
apk add sed bash make git gcc g++ gfortran tcl tcl-dev python3 python3-dev

首先把软件包管理器更新，然后用 add 安装一些依赖软件

sed 是一个文本文件修改器。会使用 vim 的用户也可以使用 vim 来更方便地修改。
bash 是终端。Alpine 默认是没有 bash 的，要自行安装。其它 Linux 都默认安装了。
git 是版本控制软件，用于拉取 OpenSees 源代码。
gcc g++ gfortran 分别是 C C++ Fortran 编译器。
tcl tcl-dev 是 tcl 依赖，如果不使用 tcl 调用 OpenSees 可以不安装。
python3 python3-dev 是 python3 开发要安装的。

编译库

下面就可以编译 OpenSees 的库了。

首先，从 GitHub 中拉取 OpenSees 的源代码，并回退到上一个发布版本。

在 Alpine 系统中，默认用户是 root 。其它操作系统一般不会以 root 登录。这里我们假设有一个用户名为 ops 。在 Alpine 中可以建立用户并切换。我们这里不建立用户了，只建立一个目录。输入

1
2
3
mkdir /home/ops
cd /home/ops
git clone https://github.com/OpenSees/OpenSees.git

现在在 ops 目录下，克隆了 OpenSees 仓库。输入 ls ，可以发现一个名为 OpenSees 的文件夹。

下面再建立两个文件夹

1
mkdir bin lib

bin 是编译 OpenSees 生成可执行文件的位置。如果不使用 tcl 这个文件夹也可以不创建。
lib 是编译成库的位置。

然后，进入 OpenSees 文件夹，回退到 v3.2.0 版本。

1
2
cd OpenSees
git checkout v3.2.0

在编译之前，要根据操作系统修改一下 Makefile.def 文件 。这里我们以 UBUNTU 为基础改动。

1
2
3
cp MAKES/Makefile.def.EC2-UBUNTU Makefile.def
sed -i 's#HOME\t\t= ./home#HOME\t\t= /home/ops#' Makefile.def
sed -i 's#/usr/lib/x86_64-linux-gnu/libtcl8.6.so#/usr/lib/libtcl8.6.so#' Makefile.def

先用 cp 语句把模板文件复制到根目录中。
sed 用于修改文本文件，可以避免交互式操作。使用 -i 可以直接在原文本中修改。这里是指把搜索到的前部分的内容更换为后面部分的内容（两部分这里用 # 分开）。

以上，编译的方式就在 Makefile 文件中定义好了。下面我们来编译即可。这里我们只把原代码编译成库，不编译二进制文件。输入

1
make libs

等待编译完成。在编译过程中会显示一些错误提示，不必理会。完成后，在 /home/ops/lib 文件夹中，可以看到生成了一系列 .a 文件。

编写分析模型

如果使用常见的 Python 调用 OpenSees 模型，速度比较慢。因此，我们采用直接写 C++ 的方法。

在 OpenSees 源码的根目录中，创建一个文件夹，命名为 mymodel 。进入，建立一个文件 main.cpp ，用于创建模型。

在 main.cpp 文件中，首先引入头文件

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#include <stdlib.h>
#include <OPS_Globals.h>
#include <StandardStream.h>
#include <ArrayOfTaggedObjects.h>
#include <Domain.h>
#include <Node.h>
#include <Truss.h>
#include <ElasticMaterial.h>
#include <SP_Constraint.h>
#include <LoadPattern.h>
#include <LinearSeries.h>
#include <NodalLoad.h>
#include <StaticAnalysis.h>
#include <AnalysisModel.h>
#include <Linear.h>
#include <PenaltyConstraintHandler.h>
#include <DOF_Numberer.h>
#include <RCM.h>
#include <LoadControl.h>
#include <BandSPDLinSOE.h>
#include <BandSPDLinLapackSolver.h>
#include <Vector.h>
#include <Python.h>

上面的头文件除 Python.h 之外，都是 OpenSees 源代码的头文件。可以在各个文件夹中找到。

那么一个模型需要哪些头文件呢？可以参考 DEVELOPER/CORE 文件夹里面的头文件。再加上自己需要使用的模型的头文件，一般就够了。

Python.h 头文件是提供 python API 的，用于使 Python 可以识别和调用，并与 Python 的数据类型相互转化。

下面建立一个函数，用于模型分析。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
Vector parametric(int fh, int a1)
{
    Domain *theDomain = new Domain();
    Node *node1 = new Node(1, 2,   0.0,  0.0);
    Node *node2 = new Node(2, 2, 144.0,  0.0);
    Node *node3 = new Node(3, 2, 168.0,  0.0);    
    Node *node4 = new Node(4, 2,  72.0, 96.0);        
    theDomain->addNode(node1);
    theDomain->addNode(node2);
    theDomain->addNode(node3);
    theDomain->addNode(node4);
    UniaxialMaterial *theMaterial = new ElasticMaterial(1, 3000);
    Truss *truss1 = new Truss(1, 2, 1, 4, *theMaterial,  a1 * 1.0);
    Truss *truss2 = new Truss(2, 2, 2, 4, *theMaterial,  5.0);    
    Truss *truss3 = new Truss(3, 2, 3, 4, *theMaterial,  5.0);        
    theDomain->addElement(truss1);
    theDomain->addElement(truss2);
    theDomain->addElement(truss3);    
    SP_Constraint *sp1 = new SP_Constraint(1, 0, 0.0, true);
    SP_Constraint *sp2 = new SP_Constraint(1, 1, 0.0, true);    
    SP_Constraint *sp3 = new SP_Constraint(2, 0, 0.0, true);
    SP_Constraint *sp4 = new SP_Constraint(2, 1, 0.0, true);    
    SP_Constraint *sp5 = new SP_Constraint(3, 0, 0.0, true);
    SP_Constraint *sp6 = new SP_Constraint(3, 1, 0.0, true);        
    theDomain->addSP_Constraint(sp1);
    theDomain->addSP_Constraint(sp2);
    theDomain->addSP_Constraint(sp3);
    theDomain->addSP_Constraint(sp4);    
    theDomain->addSP_Constraint(sp5);    
    theDomain->addSP_Constraint(sp6);    
    TimeSeries *theSeries = new LinearSeries();
    LoadPattern *theLoadPattern = new LoadPattern(1);
    theLoadPattern->setTimeSeries(theSeries);
    theDomain->addLoadPattern(theLoadPattern);
    Vector theLoadValues(2);
    theLoadValues(0) = fh * 1.0;
    theLoadValues(1) = -50.0;
    NodalLoad *theLoad = new NodalLoad(1, 4, theLoadValues);
    theDomain->addNodalLoad(theLoad, 1);
    AnalysisModel     *theModel = new AnalysisModel();
    EquiSolnAlgo      *theSolnAlgo = new Linear();
    StaticIntegrator  *theIntegrator = new LoadControl(1.0, 1, 1.0, 1.0);
    ConstraintHandler *theHandler = new PenaltyConstraintHandler(1.0e8,1.0e8);
    RCM               *theRCM = new RCM();
    DOF_Numberer      *theNumberer = new DOF_Numberer(*theRCM);    
    BandSPDLinSolver  *theSolver = new BandSPDLinLapackSolver();       
    LinearSOE         *theSOE = new BandSPDLinSOE(*theSolver);        
    StaticAnalysis    theAnalysis(*theDomain,
                  *theHandler,
                  *theNumberer,
                  *theModel,
                  *theSolnAlgo,
                  *theSOE,
                  *theIntegrator);
    theAnalysis.analyze(1);
    Vector disp = node4->getDisp();
    return disp;
}

可以看出，这段代码与 TCL 的调用方法是类似的。只不过 TCL 调用的是封装起来的方法。要想找到这些方法如何使用，只能去参考头文件。

这个函数有两个整型的变量，一个是 fh ，即水平力。一个是 a1 ，即杆 1 的面积。分析之后，输出一个 Vector 类型的 disp 。注意这里的 Vector 类型与 C++ 中的 <vector> 是不同的，它是在 OpenSees 中实现的一个类。

现在，如果再写一个 main 函数，加入对 fh 和 a1 的循环，就可以完成参数化分析了。具体方法可以参考上一篇文章 OpenSees 不同调用方法性能比较 。

现在假设我们所需要参数化分析的参数是从 Python 计算而得的。并不能直接输入给 C++ 。因此，我们把这个函数打包为一个 python 函数。

把 C++ 函数打包成 Python 库

在上文的 main.cpp 函数中，定义了一个名为 parametric 的函数，输出两个 int 型变量，输出一个 Vector 类型的值。现在我们把这个函数打包成 Python 库。

在 Python3 中，打包时应用的是 Python API ，定义于 python.h 中。打包的步骤如下

  1. 定义一个 PyObject ，为要打包的函数，把 Python 的输出类型转换为 C++ 类型，执行 C++ 函数后，再把返回的 C++ 类型转换回 Python 类型。
  2. 定义一个 PyMethodDef ，为一个数组，用于把要定义的函数名、函数文档等信息整理在一起。
  3. 定义一个 PyModuleDef ，为一个 struct 类型，用于定义 Module 的信息。
  4. 定义一个返回 PyMODINIT_FUNC 的函数，用于在 Python 引入该 Module 时初始化。

下面我们逐条来定义。

首先，定义 PyObject

1
2
3
4
5
6
7
8
9
10
static PyObject* analyze_one(PyObject* self, PyObject* args)
{
    int fh, a1;
    if (!PyArg_ParseTuple(args, "ii", &fh, &a1))
        return NULL;
    Vector v = parametric(fh, a1);
    double d1 = v(0);
    double d2 = v(1);
    return Py_BuildValue("ff", d1, d2);
}

这里我们定义一个函数，名为 analyze_one 。它在 C++ 中是一个 PyObject 的实例。在 Python 中，就是一个函数对象。

Python 把函数的变量以 Tuple 的形式传入 args 中。然后通过 PyArg_ParseTuple 这个函数来把这个 args 变量进行解析。 "ii" 是模板，表示是有两个整数的 Tuple 。

然后就可以调用 C++ 的函数，得到返回值，是 Vector 类型。这里我们把它分解为 d1 和 d2 表示两个方向的位移。

最后，使用 Py_BuildValue 把两个 double 类型的变量转化为 Python 的 float 类型，用模板 "ff" 来表示。

以上，完成了函数的封装。下面把它定义为一个 Python module 的方法。

1
2
3
4
static PyMethodDef myMethods[] = {
    {"analyze_one", analyze_one, METH_VARARGS, "Analyze one structure"},
    {NULL, NULL, 0, NULL}
};

这里我们定义一个 PyMethodDef 类型的数组 myMethods 。其中有两条记录。每条记录是一个 struct 类型。第一项为 Python 中的方法名，这里也叫 analyze_one 。第二项为上面封装的函数。第三项为一个标识，在定义一个函数时就用 METH_VARARGS 这个常量。第四项是文档。

下面再来写 module 的定义

1
2
3
4
5
6
7
static struct PyModuleDef myModule = {
    PyModuleDef_HEAD_INIT,
    "parametric",
    "Wrapped analyze module.",
    -1,
    myMethods
};

这里定义了一个 PyModuleDef 对象，是一个 struct 类型。第一项为标识，永远是用 PyModuleDef_HEAD_INIT 这个常量。第二项为 module 名，这里我们命名为 parametric 。第三项为文档，第四项一般用 -1 即可。第五项为刚才定义的方法数组。

最后，指定加载 module 时的动作即可

1
2
3
4
PyMODINIT_FUNC PyInit_parametric(void)
{
    return PyModule_Create(&myModule);
}

这里，函数名必须是 PyInit_ 加上 module 名的形式。

至此，用 Python 封装的操作就完成了。

下面我们来编译。

编译动态链接库

前文中，已经把 OpenSees 的 libs 编译完成了。下面我们把新写的 C++ 文件与它们链接起来，编译为动态链接库。

在 mymodel 文件夹中，新建一个 Makefile 文件，输入

1
2
3
4
5
6
7
include ../Makefile.def

all: parametric.so

parametric.so:
	$(LINKER) main.cpp -fPIC -shared $(INCLUDES) -L /home/ops/lib -lOpenSees \
	-I/usr/includes/python3.8  -o parametric.so

然后在输制台中，输出 make 即可编译。生成了一个 parametric.so 文件。即动态链接库文件。

使用 Python 调用

在 mymodel 目录中，建立以下 python 文件 run.py

1
2
3
4
5
6
7
8
9
10
import time
from parametric import analyze_one

startTime = time.time()
for fh in range(1, 1001):
    for a1 in range(1, 11):
        d1, d2 = analyze_one(fh, a1)
        print(d1, d2)
endTime = time.time()
print("Time:", endTime - startTime, "s")

然后，在控制台输入 python run.py 即可运行。(Alpine 系统用 python3 代替 python)

运行速度

下面来比较一下运行时间。通在我的 MacBook docker 环境中测试。测试时，为了方便，共运行3次，取3次中的中位数。由于 控制台的 I/O 输出占用时间较长，我们统计一个含输出的时间，再统计一个完全不含输出的时间 （把 print 等语句去掉）。

使用本文 C++ 封装的方法运行程序，无输出时耗时 0.561 s ，有输出时耗时 2.852 s 。相比上文的讨论，采用 OpenSeesPy 的无输出耗时 3.653 s，有输出有时 5.028 s，速度快了很多。但是与原生 C++ 的速度比（无输出 0.152 s ，有输出 0.998 s）还是要稍慢一些。

源码下载

本文和 OpenSees 不同解释器的性能比较 的源代码可以 点击下载