平面框架内力计算

function exam3_1
% 本程序为第三章的第一个算例，采用平面梁单元计算平面刚架的变形和内力
%      输入参数： 无
%      输出结果： 节点位移和单元节点力 

    PlaneframeModel ;      % 定义有限元模型
    SolveModel ;           % 求解有限元模型
    DisplayResults ;       % 显示计算结果
return ;

function PlaneframeModel
%  定义平面杆系的有限元模型
%  输入参数：
%      无
%  返回值：
%      无
%  说明：
%      该函数定义平面杆系的有限元模型数据：
%        gNode ------- 节点定义
%        gElement ---- 单元定义
%        gMaterial --- 材料定义，包括弹性模量，梁的截面积和梁的抗弯惯性矩
%        gBC1 -------- 约束条件
%        gNF --------- 集中力
%        gDF --------- 分布力

    global gNode gElement gMaterial gBC1 gNF gDF

    % 节点坐标
    %         x      y
    gNode = [0.0   0.0          % 节点 1
             0.0   4.0          % 节点 2
             3.0   0.0          % 节点 3
             3.0   4.0          % 节点 4
             4.5   4.0          % 节点 5
             6.0   0.0          % 节点 6
             6.0   4.0 ] ;      % 节点 7
     
    % 单元定义
    %          节点1  节点2  材料号
    gElement = [1      2      1       % 单元 1
                2      4      1       % 单元 2
                3      4      1       % 单元 3
                4      5      1       % 单元 4
                5      7      1       % 单元 5
                6      7      1] ;    % 单元 6
        
    % 材料性质 
    %           弹性模量  抗弯惯性矩  截面积
    gMaterial = [2.1e11   2.0e-4    1.0e-2] ;   %  材料 1

    % 第一类约束条件
    %     节点号   自由度号    约束值
    gBC1 = [ 1        1        0.0
             1        2        0.0
             1        3        0.0
             3        1        0.0
             3        2        0.0
             3        3        0.0
             6        1        0.0
             6        2        0.0
             6        3        0.0] ;

    % 集中力
    %     节点号   自由度号   集中力值
    gNF = [  5       2         -80e3] ;

    % 分布载荷（线性分布）
    %     单元号   节点1载荷值   节点2载荷值   自由度号
    gDF = [  1         -30e3         0              2
             2         -15e3       -15e3            2  ] ;
return

function SolveModel
%  求解有限元模型
%  输入参数：
%     无
%  返回值：
%     无
%  说明：
%      该函数求解有限元模型，过程如下
%        1. 计算单元刚度矩阵，集成整体刚度矩阵
%        2. 计算单元的等效节点力，集成整体节点力向量
%        3. 处理约束条件，修改整体刚度矩阵和节点力向量
%        4. 求解方程组，得到整体节点位移向量

    global gNode gElement gMaterial gBC1 gNF gDF gK gDelta

    % step1. 定义整体刚度矩阵和节点力向量
    [node_numberdummy] = size（ gNode ） ;
    gK = zeros（ node_number * 3 node_number * 3 ） ;
    f = zeros（ node_number * 3 1 ） ;

    % step2. 计算单元刚度矩阵，并集成到整体刚度矩阵中
    [element_numberdummy] = size（ gElement ） ;
    for ie=1:1:element_number
        k = StiffnessMatrix（ ie 1 ） ;
        AssembleStiffnessMatrix（ ie k ） ;
    end

    % step3. 把集中力直接集成到整体节点力向量中
    [nf_number dummy] = size（ gNF ） ;
    for inf=1:1:nf_number
        n = gNF（ inf 1 ） ;
        d = gNF（ inf 2 ） ;
        f（ （n-1）*3 + d ） = gNF（ inf 3 ） ;
    end

    % step4. 计算分布力的等效节点力，不集成到整体节点力向量中
    [df_number dummy] = size（ gDF ） ;
    for idf