MFEAOOP程序扩展：二维/三维梁单元内力计算

问题：常有小伙伴在后台私信木木，想请教在应用梁单元时，除了关心的节点位移和转角，单元内力（轴力、剪力、弯矩）也是时常关心的量，请问在程序中该怎么提取？

求解单元内力步骤：

1. 程序在求得节点位移和转角后，可根据各个单元的刚度方程求出整体坐标系下的节点力       。    
2. 此时的       是包含了各种荷载的节点力，比如重力荷载的节点等效、压强荷载的节点等效，需要减掉这些等效节点力
3. 经过以上步骤后，得到的是整体坐标系下的单元内力，也就是轴力、剪力、弯矩。需要左乘坐标旋转矩阵，才可以得到局部坐标系下的单元内力

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

相同的道理，也可以添加三维梁单元的内力计算函数。

最终的计算结果可罗列如下的效果：