题外话：因为微信的推荐机制变动，有可能大家不会第一时间看到我的文章，请大家给我的公众号标上⭐，以免错过好资源。
问题：常有小伙伴在后台私信木木，想请教在应用梁单元时，除了关心的节点位移和转角，单元内力（轴力、剪力、弯矩）也是时常关心的量，请问在程序中该怎么提取？
求解单元内力步骤：
程序在求得节点位移和转角后，可根据各个单元的刚度方程求出整体坐标系下的节点力 。 此时的 是包含了各种荷载的节点力，比如重力荷载的节点等效、压强荷载的节点等效，需要减掉这些等效节点力 经过以上步骤后，得到的是整体坐标系下的单元内力，也就是轴力、剪力、弯矩。需要左乘坐标旋转矩阵，才可以得到局部坐标系下的单元内力
MFEAOOP 是木木基于 matlab 开发的面向对象编程式的有限元通用程序，可以很方便的在原有的程序框架上进行扩充功能。比如在梁单元的内力求解上，可以在 FEASolver 类的基础上添加私有方法：
function calculateBeamForce2D(obj)
% 计算2D梁单元内力
numElements = size(obj.Element, 1);
% 初始化节点力矩阵
obj.NodeForce = zeros(numElements, 6); % 存储每个梁单元的节点力
for ie = 1:numElements
% 获取当前单元的节点编号
elementNode = obj.Element(ie, :);
elemNodeCoordinate = obj.Node(elementNode, :);
% 计算梁单元的长度 L
x = elemNodeCoordinate(:, 1);
y = elemNodeCoordinate(:, 2);
L = sqrt((x(2) - x(1))^2 + (y(2) - y(1))^2);
%% 坐标变换矩阵R
C = (x(2) - x(1)) / L;
S = (y(2) - y(1)) / L;
R = [C S 0 0 0 0;
-S C 0 0 0 0;
0 0 1 0 0 0;
0 0 0 C S 0;
0 0 0 -S C 0;
0 0 0 0 0 1];
%% 开始计算内力
% 提取当前单元的位移向量
de = zeros(6, 1);
de(1:3) = obj.U((elementNode(1) - 1) * 3 + (1:3));
de(4:6) = obj.U((elementNode(2) - 1) * 3 + (1:3));
% 使用预先存储的刚度矩阵
k = obj.ElementMatrix{ie};
% 计算局部力
localForce = k * de;
% 获取全局自由度索引
globalDOFIndex = [(elementNode(1) - 1) * 3 + (1:3), (elementNode(2) - 1) * 3 + (1:3)];
localForce = localForce - obj.FFG(globalDOFIndex);
% 计算节点力
localForce = R * localForce; % 先旋转到局部坐标系，再计算局部力
% 存储结果
obj.NodeForce(ie, :) = localForce;
end
end
相同的道理，也可以添加三维梁单元的内力计算函数。
最终的计算结果可罗列如下的效果：
参与更多互动交流，快快在下方留言区留下你的小脚印吧~
-End-
易木木响叮当
想陪你一起度过短暂且漫长的科研生活
