MFEAOOP V1.1版本正式发布啦！
自从上次将原先的 MFEA 代码修改为面向对象风格后，就一直想着添加一些好玩的功能。这不，前几天推送了几篇有关三维单元表面压强荷载如何等效为节点荷载，就想着，把这个功能移植到原有的程序中。
为了增加通用性，四面体、六面体的低阶、高阶单元都将支持表面压强的施加，这只是本次更新的一个小小一点，闲话不多说，直接先来看更新内容：
支持荷载类型：位移加载、表面压强加载（集中力我舍弃了，建议使用位移加载） 单元类型：C3D4、C3D10、C3D8、C3D20，后续的更多三维单元类型将逐步添加 内核求解：C3D20 单元增加了 14 点积分求解技术，保持精度的前提下，减少原有的高斯积分数量 稀疏矩阵存储格式：使用三元组（Triplet）形式，即行索引数组、列索引数组和非零值数组，分别对应矩阵中非零元素的位置和它们的值，大大提升了整体刚度矩阵的组装效率 线性方程组求解技术：采用直接矩阵求解技术和PCG 迭代求解技术（用户自主选择），可以大大提升刚度方程的求解效率 边界条件：使用罚函数技术计算出节点位移后,一键解出节点的支反力 前处理：支持一键提取 Abaqus-inp 文件的节点、单元、节点集、单元集、surface 等信息，供 MFEA 前处理使用 配置文件：使用 json 文件作为程序的配置文件，参数化控制材料信息、边界条件、后处理等 后处理：在 json 文件中添加后处理的控制选项，可在 Matlab 环境下进行云图绘制，可选择：场变量、单元边界颜色、面颜色、colormap、discretize 离散程度等参数 文件输出： 计算结果文件（.dat），记录各节点的坐标、位移、支反力的值 vtk 文件，可导入 Paraview 中进行可视化显示，或者使用 Python 的 Pyvista 库进行可视化也都是可以的 log 文件，记录计算各个过程的具体时间明细以及一些模型信息。
功能介绍
前处理
之前程序的前处理文件是我编制的一些数据规则供 MFEA 使用，但是从程序的通用性来讲，这样是极其不方便的，如果能继承受众面更广的商业软件的数据格式的话，用户使用起来会更顺手。
我与 Abaqus 结缘最深，那就依据拿 inp 文件开刀吧！
于是乎就手动编制了一个数据提取脚本，名字为 abaqus2MFEA.m，是不是挺通俗易懂的。程序+注释（很详尽）+警告、报错信息提示+空行总共有 940 多行...是不是看起来很多。
其实也还好，里面是由诸多功能函数组成，别看代码量大，下意识就会以为用到很多冗长的循环操作等减缓数据读取的时间，实际上这个脚本的数据解析时间非常快，昨晚测试了一个 50W 个 C3D4 单元，用这个脚本提取 inp 文件的数据，共计 5 秒不到，相比于我之前的程序，那真是快多了。
程序的具体逻辑，三两句有点说不清楚，里面还有一些小故事，可以专门写篇推文来分享。可以看一下，该脚本提取的信息有哪些：
其中，Material、ElementType、JobName、JsonName、Scale 是从 Json 文件读取的，将会在下节详细介绍。
从 inp 文件中读取的信息：Node、Element、Elset、Nset、Surface 信息，Load 和 Constr 是根据提取的信息，自动计算得出的节点边界条件信息。
Node：每一行表示节点的 x、y、z 坐标值，Element 每一行表示单元局部节点，由于是 C3D4 单元，所以是 4 列数组。
ElsetNset 数组：
包含节点集、单元集的名字、类型、具体包含哪些节点、单元。
Surface 数组：
包含：Surface 的名字、有哪些单元集（名字）以及这些单元集的面索引，在 INP 中是这样的：
*Surface, type=ELEMENT, name=Surf-1
_Surf-1_S3, S3
_Surf-1_S1, S1
_Surf-1_S4, S4
_Surf-1_S2, S2
json配置文件
JSON（JavaScript Object Notation）是一种轻量级的数据交换格式，使用键值对和数组的结构来表示数据，具有很强的可读性和可扩展性。其基本数据结构包括以下几种：
对象（Object）：由花括号
{}包围的键值对集 合，每个键值对之间用逗号分隔，键是字符串，值可以是任意的合法 JSON 类型。示例：{
"name": "Alice",
"age": 25,
"isStudent": false
}数组（Array）：由方括号
[]包围的有序值列表，每个值之间用逗号分隔，数组中的值可以是任意 JSON 类型。示例：[ "apple", "banana", "cherry" ]键值对：对象中的每个键都对应一个值。键通常是字符串类型，值可以是以下任何类型：字符串、数字、布尔值、数组、对象或
null。数据类型：
字符串（String）：用双引号括起来的文本，支持转义字符。 数字（Number）：整数或浮点数。 布尔值（Boolean）： true或false。null：表示空值或无数据。 数组（Array）：一组有序的值。 对象（Object）：键值对的无序集 合。
现将 MFEAOOP 支持的 json 语法罗列：
{
"name": "Job-1",
"scale": 1,
"inpFileName" : "../input/Abaqus/C3D20.inp",
"mesh" : {
"etype" : "C3D20"
},
"materials" : [
{
"name" : "Material-1",
"category" : "Elastic",
"data" : {
"E" : 1e6,
"v" : 0.3
}
},
{
"name" : "Material-2",
"category" : "Elastic",
"data" : {
"E" : 1e5,
"v" : 0.3
}
}
],
"bc" : [
{
"name" : "fix1",
"type": "constraint",
"direction": "xyz",
"nset" : "Set-XYZ"
},
{
"name" : "dis1",
"type": "displacement",
"direction": "xyz",
"nset" : "Set-XYZDis",
"value" : 1.0
},
{
"name" : "pressure1",
"type": "pressure",
"surface" : "Surf-1",
"value" : 1.0
},
]
"plot" : {
"fields" : ["model"],
"edgeColor" : "#000080",
"faceColor" : "#77AC30",
"colorMap" : "abaqus",
"discretize" : true,
"discretizeNum" : 12
}
}
name表示 Job 的名字，后续产生的结果文件、logw 文件、vtk 文件，都将以这个命名，类似于 Abaqus 的 Job 文件scale表示模型缩放系数inpFileName表示 inp 文件的位置，现在的前处理部分，只需要 inp 和这个配置文件即可，是不是方便许多mesh：包含单元类型，现支持 C3D4、C3D10、C3D8、C3D20，后续会不断扩充，这个模块也会根据需要来增加与网格有关的选项，定制化程度较高materials：表示材料信息，可定制多种材料，但是现在程序只支持一种材料的计算，后续如果有需要会增加材料的种类和数量bc：表示边界条件信息，支持：constraint、displacement、pressure，这三种，顾名思义也就是固定自由度、位移加载、表面压强加载。direction表示自由度方向，比如想要固定 x、y 自由度方向，就直接输入 xy，想要固定 x 方向，就直接输入 x。constraint 和 displacement都需要指定节点集！因为程序设定是直接施加在节点上的，根据 inp 文件里的节点集来就行。pressure需要指定 surfacename，也是保持和 inp 文件相同的逻辑。plot：表示绘图选项，fields为要显示的场变量信息，支持：U1、U2、U3、Umag、RF1、RF2、RF3、RFmag、model、deform，若数组内包含多个场变量，将会连着绘制多副云图，model 表示模型网格图、deform 表示变形网格图，edgeColor和faceColor表示单元边界、单元面的颜色，colorMap表示 map 信息，可以设置为 jet 等 matlab 内置的一些 colormap 或者设置为 abaqus ，即 abaqus 的云图风格，discretize表示离散程度，设置为 false 的话，就是光滑连续的，若为 true，则模型场云图呈现阶梯分布，段数由 discretizeNum 控制。
总之，按照 Abaqus 的逻辑来设计程序，以后如果有时间，完全可以将该套程序设计为一款类似 Abaqus 的 APP，让用户在使用体验上，仿佛在使用 Abaqus。
内核求解
整体刚度的存储方式采用三元组（Triplet）形式，即行索引数组、列索引数组和非零值数组，分别对应矩阵中非零元素的位置和它们的值，大大提升了整体刚度矩阵的组装效率。该处理方式，借鉴的是 Github 的开源项目: https://github.com/Qinxiaoye/FEM3D-C3D8-，将之拓展至适应所有类型的单元。
线性方程组求解时，采用直接矩阵求解技术和PCG 迭代求解技术（用户自主选择）
如此以来，在应对较大模型的求解时，求解效率也会比较理想。
边界条件的处理方式采用罚函数数值技术，在之前的推文罚函数处理边界条件（理论推导+数值实现） 里有详细分享过。该方法可以很好的处理位移加载的边界条件，可以很方便的一键求出节点的支反力，适合计算机处理。
单元测试
接下来测试几个小案例，与 Abaqus 对比一下数值精度。
模型
本次要测试的模型如下，一侧固定 XYZ 自由度方向的自由度，上表面承受表面压强荷载，以 C3D4、C3D10、C3D8、C3D20 单元为例，与 Abaqus 做对比验证 MFEAOOP程序的精度。
C3D4
C3D10
C3D8
Abaqus 的 C3D8 单元刚度矩阵做了一些修正，所以在该工况下的两者可能会出现稍微误差，后续我也会按照 Abaqus 的修正方式，将 C3D8 单元做出相应调整。注：该技术已在前期推文中对于平面等参四边形单元中数值实现。
C3D20
性能测试
在测试的案例中，网格规模都比较小，现以一个 50 万个 C3D4 单元的模型为例，看一下 MFEAOOP 计算的速度，计算日志如下：
<<< 左右滑动见更多 >>>
计算过程时间明细
| 读取 inp 文件 | 整刚组装 | 线性方程组求解 | 结果文件保存 | vtk 文件保存 | 云图绘制 | 总耗时 |
|---|---|---|---|---|---|---|
| 5.0 s | 35.5 s | 33.4 s | 3.0 s | 1.2 s | 0.2 s | 78.6 s |
读取 inp 数据的时间还是很快的，在整刚组装和线性方程组求解问题上其实还可以再进一步优化，但是高性能求解是无上限的，可能有时间为了快那十秒，需要的学习成本比想象中的大，建议有限元小白（大神除外），将主要精力放在程序的适用性、通用性、核心功能实现方面上，高性能计算可以在以上几个方面完善之后，再来精雕细琢。
