clc,clear,close all
global m b D
height=20;%桩总高
deep=9;%挖土深度
width=2;%桩间距
D=0.8;%直径
E1=3e10;%桩弹性模量
E2=3e10;%连梁弹性模量
b=0.8;%连梁高
h=0.6;%连梁宽
m=4e6;%弹簧反力系数
%%%%%%%%计算q1 q2%%%%%%%%%%%%
phi=25;%内摩擦角
q=10;%q为超载
y=19.2;%土重度
c=12;%土体粘聚力
ka=tand(45-phi/2)^2;
sigma_a=(q+y*deep)*ka-2*c*sqrt(ka);
L0=deep*tand(45-phi/2);
beta=2*width/L0-(width/L0)^2;
Pa1=beta*sigma_a;%KPa
Pa2=(1-beta)*sigma_a;%KPa
p1=-Pa1*1e3;%Pa
p2=-Pa2*1e3;%Pa
%define nodes
for i=1:height+1
    nodeCoord(i,:)=[0,i-1];%1m一个单元
end    
nodeCoord(height+2,:)=[width,height];
for i=1:height
    nodeCoord(height+2+i,:)=[width,height-i];%1m一个单元
end
node_num=length(nodeCoord(:,1));
x=nodeCoord(:,1)'; % 节点x轴方向坐标
y=nodeCoord(:,2)'; % 节点y轴方向坐标
A(1:height)=pi*D^2/4;
A(height+1)=b*h;
A(height+2:2*height+1)=pi*D^2/4;
I(1:height)=pi*D^4/36;
I(height+1)=b*h^3/12;
I(height+2:2*height+1)=pi*D^4/36;
E(1:height)=E1;
E(height+1)=E2;
E(height+2:2*height+1)=E1;
%单元信息：编号，节点1编号，节点2编号，截面面积，弹性模量
%define ele
for i=1:node_num-1
    ele(i,:)=[i,i,i+1,A(i),E(i),I(i)];
end
q1=D*p1;%桩1的线荷载
q2=D*p2;%桩2的线荷载
%载荷信息：节点编号，x方向力，y方向力，z方向力
for i=1:node_num
    if y(i)<=deep && x(i)==min(x)
        Load(i,:)=[i q1 0 0];
    elseif y(i)>deep && x(i)==min(x) 
        Load(i,:)=[i q1/(height-deep)*(height-y(i)) 0 0];
    elseif y(i)<=deep && x(i)==max(x)
        Load(i,:)=[i q2 0 0];
    elseif y(i)>deep && x(i)==max(x) 
        Load(i,:)=[i q2/(height-deep)*(height-y(i)) 0 0];
    end
end

%约束：节点编号，x方向约束，y方向约束，z方向约束
Constr=[1 NaN 0 NaN;node_num NaN 0 NaN];
Dofs=3*size(x,2); %总自由度数
EleCount=size(ele,1); %单元总数
K=zeros(Dofs,Dofs); %总体刚度矩阵
F=zeros(Dofs,1); %总体载荷列阵
U=zeros(Dofs,1); %总体位移列阵
BarLength=BarsLength(x,y,ele);
figure('Name','Undeformed Truss')
RenderFrame(ele,Load,Constr,x,y,U,1,'-k',1) %绘制桁架
figure('Name','Deformed Truss')
RenderFrame(ele,Load,Constr,x,y,U,0.5,'-k',1) %绘制变形后桁架
hold on
%遍历所有单元，将各单元刚度阵分块组装到总体刚度阵
for iEle =1:EleCount
        %该单元的两个节点的编号
        n1=ele(iEle,2);n2=ele(iEle,3);
        %计算坐标变换矩阵
        R=CoordTransform([x(n1) x(n2)],[y(n1) y(n2)],BarLength(iEle));
        %计算单元刚度矩阵 Ke=R'*ke*R;局部坐标系下的单元刚度阵转换为全局坐标下的单元刚度阵
        if y(n1)<=deep     
            ke= FrameElementKe2(ele(iEle,4),ele(iEle,5),ele(iEle,6),R,BarLength(iEle));
        else     
            ke= FrameElementKe1(ele(iEle,4),ele(iEle,5),ele(iEle,6),R,BarLength(iEle));
        end
        %将各单元刚度分块组装到总刚相应位置
        eleDof=[n1*3-2:n1*3,n2*3-2:n2*3];
        K(eleDof,eleDof)=K(eleDof,eleDof)+ke;
  
end
%形成载荷列阵--1、2、3、4自由度赋载荷值
for LoadNum=1:size(Load,1)
   for i=1:3
      F(3*Load(LoadNum)+i-3,1)=Load(LoadNum,i+1);
   end
end
%施加约束--乘大数法
for iConstr=1:size(Constr,1)
       for j=2:4
              if ~isnan(Constr(iConstr,j))
                 K(3*Constr(iConstr,1 )+j-4,3*Constr(iConstr,1 )+j-4)=...
                 1e12*K(3*Constr(iConstr,1 )+j-4,3*Constr(iConstr,1 )+j-4);
                 F(3*Constr(iConstr,1)+j-4)=1e12*Constr(iConstr,j)*...
                 K(3*Constr(iConstr,1 )+j-4,3*Constr(iConstr,1)+j-4);
              end
       end
end
U=K\F; %全局坐标系下位移,此处原来为逆但是矩阵奇异，故改成伪逆
for iEle =1:EleCount
        %计算杆局部坐标下的位移 
        n1=ele(iEle,2);n2=ele(iEle,3);
        R=CoordTransform([x(n1) x(n2)],[y(n1) y(n2)],BarLength(iEle));
        localU = R*[U(3*n1-2:3*n1,1);U(3*n2-2:3*n2,1)];
        if y(n1)<=deep
        [Ke_local] = FrameElementKeLocal1(ele(iEle,4),ele(iEle,5),ele(iEle,6),R,BarLength(iEle));
        else
        [Ke_local] = FrameElementKeLocal2(ele(iEle,4),ele(iEle,5),ele(iEle,6),R,BarLength(iEle));
        end
        EleForce(:, iEle)=Ke_local*localU; %轴力
end


MomentForce=EleForce([3,6],:);
figure
MomentForcePlot1=[-MomentForce(1,1:20)];
plot(MomentForcePlot1,[-19:0])
hold on;
MomentForcePlot2=[MomentForce(2,22:41)];
plot(-MomentForcePlot2,[0:-1:-19])
legend('前','后')
title('Moment');
figure
for i=1:node_num
    Ux(i)=U((i-1)*3+1);
end
plot(Ux(1:21),[-20:0])
hold on;
plot(Ux(22:42),[0:-1:-20])
title('Displacement');
legend('前','后')
figure
MomentForce(1,:)=-MomentForce(1,:);
ForcePlot(ele,Load,Constr,x,y,U,3,'-b',MomentForce)%绘制轴力图
title('弯矩图');