clear all;close all;clc;
%单位制mm N MPa
E=210e3;%弹性模量
t=2;%梁高度
b=100;%梁宽度
I=1/12*100*t^3;%惯性矩
L=1;%单元长度
TL=500;%梁总长度
ele_num=TL/L;%单元个数
%创建单元刚度矩阵
for i=1:ele_num
    ke(:,:,i)=BeamElementStiffness(E,I,L);%%%单元刚度矩阵的推导
end
%创建全局刚度矩阵
K=zeros((ele_num+1)*2);
for i=1:ele_num
    K=BeamAssemble(K,ke(:,:,i),i,i+1); %全局刚度矩阵
end

k=[K(3:(ele_num+1)*2-2,3:(ele_num+1)*2-2)]; %去除约束对应的刚度矩阵 
f=zeros(2*(ele_num-1),1); %节点力f
f(ele_num-1)=-400; %中间节点施加400N向下的集中里
u=k\f; %%%% inv(k)*f 高斯消元法
U=zeros((ele_num+1)*2,1);%节点位移
U(3:ele_num*2)=u;%
for i=1:ele_num
    Fe(i,:)=ke(:,:,i)*U((i-1)*2+1:(i-1)*2+4); %%%%单元力=单元刚度矩阵*节点位移
end
F=K*U; %%%%全局刚度矩阵*节点位移=节点力（所受到的不平衡反力）
figure;
x = [0:L:500]';
M=Fe(:,2);
M(ele_num+1)=-Fe(end,4);
hold on;
title('Bending Moment Diagram');
plot(x,M);
y1 = zeros(1,ele_num+1);
plot(x,y1,'k')
xlabel('length(mm)')
ylabel('Moment(N.mm)')
title('moment diagram with constant thickness');
figure;
hold on;
title('Displacment Diagram');
for i=1:length(x)
    disp_plot_constant(i)=U(2*(i-1)+1);
end
plot(x,disp_plot_constant);
save data disp_plot_constant