clear all;clc;close all;
global E u SD q0 Nx Ny
Nx=4;Ny=4;
E=2e6;%杨氏模量
u=0; %泊松比
SD=2;%二维问题
q0=-10;%每一步加载的荷载
%定义节点坐标%定义单元
Patch=1;%底部块体
[P1,E1,q0_node1,Boundary1] = Geometry(Patch, 0, 0);%返回节点单元信息、荷载施加节点、边界约束施加节点
dis_1=zeros(size(P1));
Patch=2;%上部块体
[P2,E2,q0_node2,Boundary2] = Geometry(Patch, 0, 0);
dis_2=zeros(size(P2));
%定义节点和单元
Node=[P1;P2];
Element=[E1;E2+size(P1,1)];
%绘制网格返回单元数节点数
[Nnod,Nele]=MeshPlot(Node,Element);
%定义材料刚度矩阵
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))];
time=1;
d_final=zeros(2*size(Node,1),1);

% f=figure;
% interval=5; %plot gif
while(time <100)
%     clf(f);%plot gif
    %建立下部块体刚度阵
    Patch=1;
    [K_Patch_1, F_Patch_1] = Linear_Elasticity(P1, E1, q0_node1,Boundary1,Patch);
     %建立上部块体刚度阵     
    Patch=2;    
    [K_Patch_2, F_Patch_2] = Linear_Elasticity(P2, E2,q0_node2,Boundary2, Patch);
    %将两个全局patch的刚度矩阵和荷载矩阵组装为一体
    K_Global=[K_Patch_1,zeros(size(K_Patch_1, 1),size(K_Patch_2,2)) ;zeros(size(K_Patch_2, 1),size(K_Patch_1,2)) ,K_Patch_2];
    F_Global=[F_Patch_1; F_Patch_2];       
    %考虑接触约束
    Penalty=1e7;%惩罚因子
    [K_Contact RC_force] = Contact_Formualtion(P1, P2,E1,E2, Penalty, dis_1, dis_2);
    K_Global=K_Global+K_Contact;
    d=linsolve(K_Global, F_Global);
    d1=d(1:SD*size(P1,1));
    d2=d(SD*size(P1,1)+1:SD*(size(P1,1)+size(P2,1)));
    dis_1=dis_1+transpose(reshape(d1, size(dis_1,2), size(dis_1,1)));
    dis_2=dis_2+transpose(reshape(d2, size(dis_2,2), size(dis_2,1)));
   	d_final=d_final+d;  

    time=time+1; 
    
    
    
%     %plot gif 绘制GIF动图
%     [NodeStrain,NodeStress]=CalculateStrainAndStress2(d_final,D,Node,Element);
%     PlotContour(Node,Element,d_final,NodeStress(2,:))%sigma_11   
%     axis([0 1 -0.1 0.7])
%     caxis([-1591 123])
%     colorbar off
%     drawnow;
%     F=getframe(f);
%     I=frame2im(F);
%     [I,map]=rgb2ind(I,256);
%     if time == 2
%         imwrite(I,map,'sigma_22.gif','gif', 'Loopcount',inf,'DelayTime',0);
%     elseif mod(time,interval)==0
%         imwrite(I,map,'sigma_22.gif','gif','WriteMode','append','DelayTime',0);
%     end
  
end
% 存储结点位移
gNTu = [dis_1;dis_2];

% 绘制变形图
SF=1;
PlotResults2(Node,Element,gNTu,100);
[NodeStrain,NodeStress]=CalculateStrainAndStress2(d_final,D,Node,Element);
figure
PlotContour(Node,Element,d_final,NodeStress(1,:))%sigma_11
title('\sigma11');
figure
PlotContour(Node,Element,d_final,NodeStress(2,:))%sigma_22
title('\sigma22');