%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [P,E,N_ID_p1,Boundary] = Geometry(Patch, dis_1, dis_2)
global  Nx Ny
    global time;
    
    gap=-0.0001;%两物体间隙

    switch Patch %options of Patch
        case 1 %下部块体                 
            P=[0,0; 0.25,0; 0.5,0;  0.75,0; 1,0;  
               0,0.125; 0.25,0.125; 0.5,0.125;  0.75,0.125; 1,0.125;
               0,0.25; 0.25,0.25; 0.5,0.25;  0.75,0.25; 1,0.25;
               0,0.375; 0.25,0.375; 0.5,0.375;  0.75,0.375; 1,0.375;
               0,0.5; 0.25,0.5; 0.5,.50;  0.75,0.5; 1,0.5
               ];%节点坐标
           
           P=P+ dis_1;%通过位移更新节点坐标
           %定义单元
           ele_num=0;
           for i=1:Nx
               for j=1:Ny
                   ele_num=ele_num+1;
                   E(ele_num,:)=[(i-1)*(Nx+1)+1+(j-1),(i-1)*(Nx+1)+2+(j-1),i*(Nx+1)+2+(j-1),i*(Nx+1)+1+(j-1)];
               end
           end
           %约束边界条件
           %边界条件：type = 1 -- fixed disp [ 1--节点固定; 0--节点自由 ]; 
            %          type = 2 -- point load 
            % node, x, y, z, rx, ry, rz, type 
            for i=1:Nx+1
                Boundary(i,:) =  [i 1.0 1.0 1.0 0.0 0.0 0.0 1];
            end
            N_ID_p1=[];
        case 2%上部块体
        
            P=[0.2,0.5; 0.3,0.5; 0.4,0.5;  0.5,0.5; 0.6,0.5;  
               0.2,0.55; 0.3,0.55; 0.4,0.55;  0.5,0.55; 0.6,0.55;
               0.2,0.6; 0.3,0.6; 0.4,0.6;  0.5,0.6; 0.6,0.6;
               0.2,0.65; 0.3,0.65; 0.4,0.65;  0.5,0.65; 0.6,0.65;
               0.2,0.7; 0.3,0.7; 0.4,0.7;  0.5,0.7; 0.6,0.7
               ];
            P(:,2)=P(:,2)+gap; %create contact gap between two patches
            P(:,1)=P(:,1);
            %update the coordinates using displacement         
           P=P+dis_2;
            %定义单元
           ele_num=0;Nx=4;Ny=4;
           for i=1:Nx
               for j=1:Ny
                   ele_num=ele_num+1;
                   E(ele_num,:)=[(i-1)*(Nx+1)+1+(j-1),(i-1)*(Nx+1)+2+(j-1),i*(Nx+1)+2+(j-1),i*(Nx+1)+1+(j-1)];
               end
           end
           %定义分布荷载的单元节点编号
            E_ID_p1 = Nx*(Ny-1)+1:Nx*Ny;      % 分布荷载作用单元编号
            N_ID_p1 = zeros(size(E_ID_p1,2),size(E,2));%分布荷载所在单元节点编号
            j = 1; 
            for i_q = E_ID_p1
                 N_ID_p1(j,:) = E(i_q,:);              %分布荷载所在单元节点编号
                 j = j+1;
            end
           %约束边界条件
           %边界条件：type = 1 -- fixed disp [ 1--节点固定; 0--节点自由 ]; 
            %          type = 2 -- point load 
            % node, x, y, z, rx, ry, rz, type 
            for i=1:Nx+1
                Boundary(i,:) =  [(Nx+1)*Ny+i 1.0 0.0 0.0 0.0 0.0 0.0 1];%保证上部块体沿着竖向运动，防止矩阵奇异
            end

                                      
    end

end