MFEAOOP V1.4版本震撼来袭！更丰富的单元库！
1.4版本更新内容
增加单元类型： 轴对称单元: CAX3、CAX4、CAX6、CAX8、CAX8R，均支持表面压强荷载和重力荷载 考虑剪切变形的 Euler-Bernoulli 平面梁单元B23M 增加89种colormap配色方案：matlab内置+Matplotlib官网+自定义配色方案 增加梁单元内力(轴力、剪力、弯矩)计算输出 二维梁：Fx Fy RM3 三维梁：Fx Fy Fz RM1 RM2 RM3 所有单元类型，除了梁单元，其余均采用mex文件加速，组装整体刚度矩阵的时间可以提升一倍多，梁单元因为截面类型不同，所需参数不同，不利于程序的统一实现，因此不采用mex文件加速 后处理部分增加节点标号显示： "nodeLabel" : true增加MPC多点约束功能（对应于Abaqus的TIE功能，即指定节点的自由度被约束）
MFEAOOP V1.4所有内容
支持荷载类型： 位移荷载 表面压强荷载 重力荷载 单元类型（共计支持39种单元类型）： C3D4、C3D10、C3D8、C3D20、C3D8R、C3D20R、C3D8I、C3D8Bar、C3D8SR CAX3、CAX4、CAX6、CAX8、CAX8R CPS4、CPS8、CPS3、CPS6、CPS4R、CPS8R、CPS4I、CPS4Bar、CPS4SR CPE4、CPE8、CPE3、CPE6、CPE4R、CPE8R、CPE4I、CPSEBar、CPE4SR T2D2、T3D2 B21、B23、B31、B33、B23M 杆系单元： 平面单元： 轴对称单元： 空间实体单元： 整体刚度矩阵组装：稀疏矩阵存储格式，使用三元组（Triplet）形式，即行索引数组、列索引数组和非零值数组。使用mex文件加速，组装整体刚度矩阵的时间可以提升一倍多。 线性方程组求解技术：提供直接求解法和pcg迭代求解器选项，面对较大模型时可以大大提升刚度方程的求解效率 边界条件：使用罚函数技术 节点应力磨平处理: 对于单元形状不均匀的结构模型，对于共节点处的应力采用体积/面积加权磨平处理 对于单元形状均匀的结构模型，对于共节点处的应力采用直接绕节点平均磨平处理 支持多种材料属性的赋予 场输出:U、RF、UR、RFM、S 前处理：解析 Abaqus-inp 文件： Node节点坐标 Element单元局部节点编号 Nset节点集 Elset单元集 Surface单元面信息 配置文件：使用 json 文件作为程序的配置文件 材料信息 截面属性 边界条件 求解器 绘图控制 文件输出： 结果文件（.dat），记录各节点位移、支反力等 vtk 文件，可导入 Paraview 中进行可视化显示，或者使用 Python 的 Pyvista 库进行可视化也都是可以的 log 文件，记录计算各个过程的具体时间明细以及一些模型信息。
MPC多点约束
MPC约束节点集的格式如下：
{
"name": "fix2",
"type": "mpc",
"couplings": [
{
"master_nodes": "Set-M", // 约束节点集 合1
"slave_nodes": "Set-S", // 约束节点集 合2
"dof": [1, 2] // 约束的自由度（说明1自由度和2自由度均被约束）
},
{
"master_nodes": "Set-p11",
"slave_nodes": "Set-p12",
"dof": [1]
}
]
}
上面的约束条件表示：约束节点集 合Set-M和节点集 合Set-S的1、2自由度进行TIE连接（变形过过程中节点集Set-M和Set-S的1、2自由度值保持一致），约束节点集 合Set-p11和节点集 合Set-p12的1自由度进行TIE连接。
接下来以一个平面模型为例，展示程序的MPC功能。
程序中将11-12号节点以及4-5号节点的1、2自由度进行TIE连接，位移云图中，这两对节点值是完全一样的。
mex文件加速
在MATLAB中，MEX文件是一种利用C、C++或Fortran编写的外部函数，可以在MATLAB中直接调用。这种文件具有显著的加速效果，尤其适合需要大量计算或重复运算的程序。通过MEX文件，将计算密集型的代码段编译为机器码，MATLAB可以直接调用这些代码，而无需解释执行，从而获得接近原生C/C++程序的执行速度。
将程序中的所有单元刚度矩阵计算函数一起转换为mex文件，可编制以下脚本：
clc; clear
% 梁单元不进行mex处理，因为截面参数多变，不利于统一封装
% 设置基本配置
cfg = coder.config('mex');
cfg.EnableOpenMP = true; % 启用多线程支持
% 定义单元类型列表和相应的节点数量
elementTypes = {
'T2D2', 'T3D2', ...
'CPS4', 'CPS8', 'CPS3', 'CPS6', 'CPS4R', 'CPS8R', 'CPS4I', 'CPS4Bar', 'CPS4SR', ...
'CPE4', 'CPE8', 'CPE3', 'CPE6', 'CPE4R', 'CPE8R', 'CPE4I', 'CPE4Bar', 'CPE4SR', ...
'C3D4', 'C3D10', 'C3D8', 'C3D20', 'C3D8R', 'C3D20R', 'C3D8I', 'C3D8Bar', 'C3D8SR'
};
% 定义默认输入参数结构体示例
% 示例材料属性
% 遍历每种单元类型并生成 MEX 文件
for i = 1:numel(elementTypes)
elementType = elementTypes{i};
% 根据不同的单元类型设置节点数和维度
switch elementType
case {'T2D2'}
materialProps = struct('E', 210e9, 'A', 0.3);
nodeCoordinates = zeros(2, 2); % 2D 桁架单元，2 个节点，每节点 2 个坐标
case {'T3D2'}
materialProps = struct('E', 210e9, 'A', 0.3);
nodeCoordinates = zeros(2, 3); % 3D 桁架单元，2 个节点，每节点 3 个坐标
case {'CPS4', 'CPS4R', 'CPS4I', 'CPS4Bar', 'CPS4SR', ...
'CPE4', 'CPE4R', 'CPE4I', 'CPE4Bar', 'CPE4SR'}
materialProps = struct('E', 210e9,'v', 0.3,'t', 0.3);
nodeCoordinates = zeros(4, 2); % 2D 四节点平面单元，4 个节点，每节点 2 个坐标
case {'CPS8', 'CPS8R'}
materialProps = struct('E', 210e9,'v', 0.3,'t', 0.3);
nodeCoordinates = zeros(8, 2); % 2D 八节点平面单元，8 个节点，每节点 2 个坐标
case {'CPS3', 'CPE3'}
materialProps = struct('E', 210e9,'v', 0.3,'t', 0.3);
nodeCoordinates = zeros(3, 2); % 2D 三节点平面单元，3 个节点，每节点 2 个坐标
case {'CPS6', 'CPE6'}
materialProps = struct('E', 210e9,'v', 0.3,'t', 0.3);
nodeCoordinates = zeros(6, 2); % 2D 六节点平面单元，6 个节点，每节点 2 个坐标
case {'CPE8', 'CPE8R'}
materialProps = struct('E', 210e9,'v', 0.3,'t', 0.3);
nodeCoordinates = zeros(8, 2); % 2D 八节点平面单元，8 个节点，每节点 2 个坐标
case 'C3D4'
materialProps = struct('E', 210e9, 'v', 0.3);
nodeCoordinates = zeros(4, 3); % 3D 四节点四面体单元，4 个节点，每节点 3 个坐标
case 'C3D10'
materialProps = struct('E', 210e9, 'v', 0.3);
nodeCoordinates = zeros(10, 3); % 3D 十节点四面体单元，10 个节点，每节点 3 个坐标
case {'C3D8', 'C3D8R', 'C3D8I', 'C3D8Bar', 'C3D8SR'}
materialProps = struct('E', 210e9, 'v', 0.3);
nodeCoordinates = zeros(8, 3); % 3D 八节点六面体单元，8 个节点，每节点 3 个坐标
case {'C3D20', 'C3D20R'}
materialProps = struct('E', 210e9, 'v', 0.3);
nodeCoordinates = zeros(20, 3); % 3D 二十节点六面体单元，20 个节点，每节点 3 个坐标
otherwise
warning('Unknown element type: %s. Skipping...', elementType);
continue;
end
try
eval(['codegen -config cfg ' elementType ' -args {materialProps, nodeCoordinates};']);
fprintf('Successfully generated MEX file for %s.\n', elementType);
catch ME
fprintf('Failed to generate MEX file for %s: %s\n', elementType, ME.message);
end
end
从以上程序可以看到在转换过程中，需要明确给定形参的初始值，在设计刚度矩阵时，我均采用这样的形参，以CPS4单元为例：
function k=CPS4(prop, elemNodeCoordinate)
DOF = 2;
elenode = length(elemNodeCoordinate);
k = zeros(DOF*elenode,DOF*elenode);
t = prop.t;
stressType = 'planeStress';
D = Dmatrix2D(prop, stressType);
ngp = 4;
[weights,locations]=gauss2D(ngp);
for n=1:size(weights,1)
xi_Gauss=locations(n,1);
eta_Gauss=locations(n,2);
[B, Jacob] = Bmatrix2D(elenode, elemNodeCoordinate, DOF, xi_Gauss, eta_Gauss);
k = k+B.'*D*t*B*weights(n)*det(Jacob);
end
end
其中prop是从json读取的材料参数，是个可变参数，单元类型不同，用到的参数也不同，平面单元中用到的是弹性模量、泊松比、厚度，所以我在设计：materialProps = struct('E', 210e9,'v', 0.3,'t', 0.3)。
轴对称单元
为丰富程序的单元库，本次更新特意新添加了五种轴对称单元：CAX3、CAX4、CAX6、CAX8、CAX8R，具体的理论部分我会找时间更新至《有限元基础编程百科全书》，以CAX8单元为例，工况如下：
大家可能注意到MFEAOOP的应力场有些怪怪的，即x=0的轴上为空白，这是为什么呢？程序刚运行出这样的结果的时候我也会很惊讶，后来仔细思考，发现这是合理的，因为我在节点应力计算时，我直接将节点的坐标代入应力计算公式，而不是高斯积分点，所以会导致x=0处的节点应力出现奇异的情况。
这个时候可以更改应力计算方案，将高斯积分点代入应力计算公式，然后进行外插即可，但是这样以来，程序需要根据不同的单元类型进行选择合适的形函数，多少有些繁琐。
我采用的方法时，将x=0的节点x坐标设为0.0001，做一个小小的修正，即可以避免数值奇异的情况，效果如下：
colormap
为增加后处理的灵活性，我将matplotlib的colormap添加到MFEAOOP中，供用户使用，效果如下：
bug修复
梁单元在约束自由度时，比如要约束全部自由度，需要xyzUR1UR2UR3
程序心得
编码格式要统一 尽量使用英文报错、警告信息 程序出现bug，不要慌，根据报错信息，查询错误源头，如果不会高级调试就使用最简单的方式：打印变量 如果不想再看程序了，那就关闭程序，隔几天后再来编写你的程序，这样你的思路会清晰很多 不追求华丽花哨的代码，追求简洁、易读、易维护的代码
由于添加的单元类型较多，精度验证需要占据大量篇幅，没有贴在本次推文中，不过大可放心，每种单元类型在添加时均和Abaqus做了验证分析，均验证通过！感兴趣的小伙伴可逐个验证。
程序目前先不继续更新功能了，接下来就是补充完善理论部分，及时更新《有限元基础编程百科全书》，程序的思维导图也在绘制中，不久将全部公布，感谢大家关注。
有关程序更新的功能就介绍到这里，感谢您的阅读。整套程序已发布在知识星球中，后台回复：星球，即可加入，对源程序的疑问，可在星球内详细讨论。
