clear all;close all;clc;
EI=1;
L=1;%杆件长度
Ele_L=0.01;%单元尺寸
ele_num=L/Ele_L;%单元个数

for i=1:ele_num
    ke(:,:,i)=BeamElementStiffness(EI,Ele_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,3:(ele_num+1)*2)]; %***去除约束对应的刚度矩阵 

f=zeros(2*(ele_num),1); %节点力f

for i=1:ele_num
    f(2*(i-1)+1)=1*Ele_L;%****分布力大小为1*单元长度就是每个单元节点上的节点力，均布荷载转节点力
end


u=k\f;
U=zeros((ele_num+1)*2,1);%节点位移
U(3:(ele_num+1)*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 : Ele_L:L]';
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')
ylabel('Moment')
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);
figure;
hold on;
title('Rotation Diagram');
for i=1:length(x)
    rot_plot_constant(i)=U(2*(i));
end
plot(x,rot_plot_constant);
save data disp_plot_constant