% 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];
% fv=[32,1 100;33 2 -50]
% fv=[32,1 400000]
% fv=[33,2 -33333.3333]
% t=20;
% E=2e5;
% Nu=0.25;
flagfilename=1;
while flagfilename > 0
inputfilename=uigetfile('*.xls;*.xlsx','Select xls file as FEM Input');
if inputfilename > 0, flagfilename=0;, end
end

ncoor=xlsread(inputfilename,'ncoor');
nconn=xlsread(inputfilename,'nconn');
t=xlsread(inputfilename,'thickness');


E=xlsread(inputfilename,'material','A1');
Nu=xlsread(inputfilename,'material','B1');



nel=size(nconn,1)

nnodes=size(ncoor,1);



KK=zeros(2*nnodes,2*nnodes);

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);
flagfv=1;
flagbc=1;
flagscale=1;


while flagfv+flagbc+flagscale > 0
    
    if flagscale==1,
        deformscaleX=xlsread(inputfilename,'scale','A1');
        deformscaleY=xlsread(inputfilename,'scale','B1');
    end
    
    if flagfv==1 || flagbc==1,
        if flagfv==1,
            fv=xlsread(inputfilename,'force');
            nfv=size(fv,1);
            F((fv(:,1)-1)*2+fv(:,2))=fv(:,3) 
        end
        if flagbc==1,
            bc=xlsread(inputfilename,'bc');
            nbc=size(bc,1);
            index_bcdof=[(bc(:,1)-1)*2+bc(:,2)]
            index_ugdof=[1:2*nnodes]';
            index_rugdof=setdiff(index_ugdof,index_bcdof)
        end
    
    
        rKK=KK(index_rugdof, index_rugdof)
        rF=F(index_rugdof);
    
        rU=inv(rKK)*rF;
        U=zeros(nnodes*2,1);
        U(index_rugdof)=rU
        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);
            qptsxy(i,:,:)=qpxy;
            sss(i,:,:)=s';
        end        
    end
    
    figure(1)
    clf;
    title('Deformed shape and indication Ux');
    xdata=reshape( (ncoor(nconn(:,:),1)+deformscaleX*U(index_ugdof(2*(nconn(:,:)-1)+1,1))), size(nconn,1),size(nconn,2));
    ydata=reshape( (ncoor(nconn(:,:),2)+deformscaleX*U(index_ugdof(2*(nconn(:,:)-1)+2,1))), size(nconn,1),size(nconn,2));
    udata=reshape( U(index_ugdof(2*(nconn(:,:)-1)+1,1)), size(nconn,1),size(nconn,2));
    patch(xdata',ydata',udata');
    colorbar;
    
    
    figure(2)
    clf;
    title('Deformed shape and indication Uy');
    vdata=reshape( U(index_ugdof(2*(nconn(:,:)-1)+2,1)), size(nconn,1),size(nconn,2));
    patch(xdata',ydata',vdata');
    
    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))]);
    title('Normal Stress - Sigma x @ Gaussian Points');
    
    figure(5)
    clf;
    colorbar;
    axis([min(ncoor(:,1)) max(ncoor(:,1)) min(ncoor(:,2)) max(ncoor(:,2))]);
    title('Normal Stress - Sigma y @ Gaussian Points');
    
    figure(6)
    clf;
    colorbar;
    axis([min(ncoor(:,1)) max(ncoor(:,1)) min(ncoor(:,2)) max(ncoor(:,2))]);
    title('Shear Stress - Tau xy @ Gaussian Points');
    

        
    figure(4)
    patch(qptsxy(:,:,1)', qptsxy(:,:,2)', sss(:,:,1)');
    figure(5)
    patch(qptsxy(:,:,1)', qptsxy(:,:,2)', sss(:,:,2)');
    figure(6)
    patch(qptsxy(:,:,1)', qptsxy(:,:,2)', sss(:,:,3)');
    choice = menu('Choose a choice to modify','scale','Force','Boundary conditions','scale + Force', 'scale + Boundary conditions', 'Force + Boundary conditions','scale + Force + Boundary conditions','CLOSE')
    flagscale=0;
    flagfv=0;
    flagbc=0;
    switch choice
        case {0, 8}
            flagscale=0;
            flagfv=0;
            flagbc=0;
        case 1,
            flagscale=1;
        case 2,
            flagfv=1;
        case 3,
            flagbc=1;
        case 4,
            flagscale=1;
            flagfv=1;
        case 5,
            flagscale=1;
            flagbc=1;
        case 6,
            flagfv=1;
            flagbc=1;
        case 7,
            flagscale=1;
            flagfv=1;
            flagbc=1;
        otherwise,
            flagscale=0;
            flagfv=0;
            flagbc=0;
    end
    if flagfv+flagbc+flagscale > 0
        flagfilename=1;
        while flagfilename > 0
            inputfilename=uigetfile('*.xls;*.xlsx','Select xls file as FEM Input');
            if inputfilename > 0, flagfilename=0;, end
        end
    end

end        
disp('GOOD BYE');