clear;clc;;close all
format short
E=210e9;%弹性模量，单位Pa
u=0.4;%泊松比，无量纲
q=100;%外载荷（面力），单位N/㎡
t=0.1;%厚度，单位m
D=[E/(1-u^2),u*E/(1-u^2),0;u*E/(1-u^2),E/(1-u^2),0;0,0,E/(2*(1+u))];%应力应变矩阵
bt=(3+0.5407762)*10^11;%乘大数法的大数
long=5;%长5m
height=1;%高1m
ele_size=[0.2,0.2];%单元的两个直角边尺寸；5m长，十个单元,1m高，4个单元
ele_H=long/ele_size(1);%水平单元数
ele_V=height/ele_size(2);%竖向单元数
ele_num=ele_H*ele_V*2;%单元数
node_num=(ele_H+1)*(ele_V+1);%节点数
DOF=2;%每个节点自由度

%定义每个单元的节点号
num=0;
for i=1:ele_H
    for j=1:ele_V
        num=num+1;
        element(num,1)=1+(j-1)+(i-1)*(ele_V+1);
        element(num,2)=1+(j-1)+i*(ele_V+1);
        element(num,3)=2+(j-1)+(i-1)*(ele_V+1);
        num=num+1;
        element(num,1)=1+j+(i-1)*(ele_V+1);
        element(num,2)=1+(j-1)+i*(ele_V+1);
        element(num,3)=2+(j-1)+i*(ele_V+1);;
    end
end

%定义节点坐标
num=0;
for i=1:ele_H+1
    for j=1:ele_V+1
        num=num+1;
        node(num,1)=(i-1)*ele_size(1);%x
        node(num,2)=(j-1)*ele_size(2);%y
    end
end
        
A=1/2*ele_size(1)*ele_size(2);%单元面积5/10*1/4矩形
%每个单元自由度在整体矩阵中的索引位置
ele_DOF=zeros(ele_num,DOF*3);
for i=1:ele_num
    for j=1:length(element(i,:))%三角形三个节点循环
        ele_DOF(i,1+2*(j-1):2+2*(j-1))=2*(element(i,j)-1)+1:2*(element(i,j)-1)+2;
    end
end

KZ=zeros(node_num*2,node_num*2);
for i=1:ele_num 
    ke=ele_stiff_matrix(element,node,i,A,D,t);%单元刚度矩阵
    KZ(ele_DOF(i,:),ele_DOF(i,:)) = KZ(ele_DOF(i,:),ele_DOF(i,:))+ke;
end

P=[]; P(2*(ele_V+1))=-1/2*q*t*long/ele_H;
for i=1:ele_H
      P(10+2*(ele_V+1)*(i-1))=-q*t*long/ele_H;%5/10是每个单元的长度
end
P(node_num*2)=-1/2*q*t*5/10;
PZ=P';%总结点载荷列阵
%乘大数法施加约束条件，最左侧节点全部固定
cons_DOF=[1:(ele_V+1)*2];%约束自由度编号
Disp=zeros(length(cons_DOF),1);%位移为0
for i=1:length(cons_DOF)
      KZ(cons_DOF(i),cons_DOF(i))=bt*KZ(cons_DOF(i),cons_DOF(i));
      PZ(cons_DOF(i))=Disp(i)*KZ(cons_DOF(i),cons_DOF(i));
end

a=inv(KZ)*PZ;%结点位移列阵
for i=1:ele_num%80个单元
      sgm(:,i)=stress_calculate(element,node,i,A,D,t,a);%单元应力
end
  
%重新排列位移矩阵
SF=1/5*1/max(abs(a));%位移放大系数
new_node=node+SF*reshape(a,2,length(a)/2)';%计算变形后的坐标

figure
patch('Faces',element,'Vertices',node,'FaceColor','none','EdgeColor','r')
hold on;
patch('Faces',element,'Vertices',new_node,'FaceColor','none')
xlabel('x-coordinate')
ylabel('y-coordinate')
title('deformation')
axis equal
%应力  
figure
patch('Faces',element,'Vertices',node,'FaceVertexCData',sgm(1,:)','FaceColor','flat')
axis equal
xlabel('x-coordinate')
ylabel('\sigma (Pa)')
title('\sigma_{11}')
colorbar
colormap(jet)

% figure
% patch('Faces',element,'Vertices',node,'FaceVertexCData',sgm(2,:)','FaceColor','flat')
% axis equal
% colorbar
% colormap(jet)
% title('sigma22')

% figure
% patch('Faces',element,'Vertices',node,'FaceVertexCData',sgm(3,:)','FaceColor','flat')
% xlabel('x-coordinate')
% ylabel('\sigma (Pa)')
% axis equal
% colorbar
% colormap(jet)
% title('sigma12')