deformscaleX=100;
deformscaleY=100;
ncoor=[0 0; 0 100;0 200;40 0;40 100;40 200;80 0;80 100; 80 200;120 0;120 100;120 200;160 0;160 100;160 200;200 0; 200 100; 200 200;240 0;240 100;240 200;280 0;280 100;280 200;320 0;320 100;320 200;360 0;360 100;360 200;400 0;400 100; 400 200];
nconn=[1 4 5 2;4 7 8 5;7 10 11 8;10 13 14 11; 13 16 17 14; 16 19 20 17; 19 22 23 20; 22 25 26 23;25 28 29 26;28 31 32 29;2 5 6 3;5 8 9 6;8 11 12 9; 11 14 15 12; 14 17 18 15;17 20 21 18;20 23 24 21; 23 26 27 24; 26 29 30 27; 29 32 33 30];
nnodes=size(ncoor,1);
nel=size(nconn,1)

bc=[1 1;1 2;2 1;2 2; 3 1;3 2];
nbc=size(bc,1);

%fv=[32,1 100;33 2 -50]
%fv=[32,1 400000]
fv=[33,2 -33333.3333]
nfv=size(fv,1);


KK=zeros(66,66);
t=20;
E=2e5;
Nu=0.25;
for i=1:nel
    disp(i)
    X= ncoor(nconn(i,:),1)';
    Y= ncoor(nconn(i,:),2)';
    k=Q4_2Dk(X,Y,t,E,Nu);
    s1=[2*nconn(i,:)-1; 2*nconn(i,:)]
    sgdof=s1(:);
    
    KK(sgdof, sgdof)=KK(sgdof,sgdof) + k;
end
    F=zeros(nnodes*2,1);
    
    F((fv(:,1)-1)*2+fv(:,2))=fv(:,3)    
    
    index_bcdof=[(bc(:,1)-1)*2+bc(:,2)]
    index_ugdof=[1:66]';
    index_rugdof=setdiff(index_ugdof,index_bcdof)
    rKK=KK(index_rugdof, index_rugdof)
    rF=F(index_rugdof)
    
    rU=inv(rKK)*rF;
    U=zeros(nnodes*2,1);
    U(index_rugdof)=rU
    
    figure(1)
    clf;
    for i=1:nel
        disp(i)
        patch( (ncoor(nconn(i,:),1)+deformscaleX*U(index_ugdof(2*(nconn(i,:)-1)+1,1)))', (ncoor(nconn(i,:),2)+deformscaleX*U(index_ugdof(2*(nconn(i,:)-1)+2,1)))', U(index_ugdof(2*(nconn(i,:)-1)+1,1)))
    end
    colorbar;
    figure(2)
    clf;
    for i=1:nel
        disp(i)
        patch( (ncoor(nconn(i,:),1)+deformscaleX*U(index_ugdof(2*(nconn(i,:)-1)+1,1)))', (ncoor(nconn(i,:),2)+deformscaleX*U(index_ugdof(2*(nconn(i,:)-1)+2,1)))', U(index_ugdof(2*(nconn(i,:)-1)+2,1)))
    end
    colorbar;
%    figure(3)
%    clf;
%    quiver(ncoor(:,1),ncoor(:,2),U(1:2:66),U(2:2:66))
%    pause(4)
    
    figure(4)
    clf;
    colorbar;
    axis([min(ncoor(:,1)) max(ncoor(:,1)) min(ncoor(:,2)) max(ncoor(:,2))])
    
    figure(5)
    clf;
    colorbar;
    axis([min(ncoor(:,1)) max(ncoor(:,1)) min(ncoor(:,2)) max(ncoor(:,2))])
    
    figure(6)
    clf;
    colorbar;
    axis([min(ncoor(:,1)) max(ncoor(:,1)) min(ncoor(:,2)) max(ncoor(:,2))])
    
    for i=1:nel
        disp(i)
        X= ncoor(nconn(i,:),1)';
        Y= ncoor(nconn(i,:),2)';
        s1=[2*nconn(i,:)-1; 2*nconn(i,:)]
        sgdof=s1(:);
        [s,e,qpxy]=Q4_2DStress_Strain(X,Y,U(sgdof),E,Nu);
        figure(4)
        patch(qpxy(:,1)', qpxy(:,2)', s(1,:));
        figure(5)
        patch(qpxy(:,1)', qpxy(:,2)', s(2,:));
        figure(6)
        patch(qpxy(:,1)', qpxy(:,2)', s(3,:));
    end
    