clear;clc;;close all
tic
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^51;%乘大数法的大数
DOF=2;%每个节点自由度

%定义每个单元的节点号
fname='M1.mphtxt';
[node, element] = Readmesh( fname );
%定义每个单元的节点号
for i=1:length(element(:,1))
    x1=node(element(i,1),1);
    x2=node(element(i,2),1);
    x3=node(element(i,3),1);
    y1=node(element(i,1),2);
    y2=node(element(i,2),2);
    y3=node(element(i,3),2);
    A(i)=(x1*y2+x2*y3+x3*y1-x1*y3-x2*y1-x3*y2);%单元面积S=(x1y2-x1y3+x2y3-x2y1+x3y1-x2y2)
end

ele_num=length(element(:,1));%单元数
node_num=length(node(:,1));%节点数
%每个单元自由度在整体矩阵中的索引位置
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

%find constrained side and force side
num=0;
num2=0;
num3=0;
for i=1:node_num
    min_x=min(node(:,1));
    max_x=max(node(:,1));
    min_y=min(node(:,2));
    if node(i,1)==min_x
        num=num+1;
        const_side_y(num)=i;%左边界
    elseif node(i,1)==max_x
        num2=num2+1;
        force_side(num2)=i;%右边界
    elseif node(i,2)==min_y
        num3=num3+1;
        const_side_x(num3)=i;%下边界
    end
end     
        
    
P=zeros(1,node_num*2);
for i=1:length(force_side)
      P(force_side(i)*2-1)=q;%5/10是每个单元的长度
end
PZ=P';%总结点载荷列阵
%乘大数法施加约束条件，最左侧节点全部固定
for i=1:length(const_side_y)
    cons_DOF_y(i)=[const_side_y(i)*2-1];%左边界对称边界条件
end
for i=1:length(const_side_x)
    cons_DOF_x(i)=[const_side_x(i)*2];%下边界对称边界条件
end

cons_DOF=[cons_DOF_y,cons_DOF_x];
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

KZ_xishu=sparse(KZ);%为了加快计算速度，采用稀疏矩阵进行计算
% a=KZ_xishu\PZ;%结点位移列阵
a=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','EdgeColor','none')
axis equal
xlabel('x-coordinate')
ylabel('\sigma (Pa)')
title('\sigma_{11}')
colorbar
colormap(jet)
%验证圆孔出1：3的mises应力比
node_ID1=10818; %(0.1,0)
node_ID2=12373;%(0,0.1)
ele1=find(element(:,2)==node_ID1);
ele2=find(element(:,2)==node_ID2);
mises1=sqrt((sgm(1,ele1)+sgm(2,ele1))^2-3*(sgm(1,ele1)*sgm(2,ele1)-sgm(3,ele1)^2));
mises2=sqrt((sgm(1,ele2)+sgm(2,ele2))^2-3*(sgm(1,ele2)*sgm(2,ele2)-sgm(3,ele2)^2));
ratio=mises2/mises1
toc
% 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')