% 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');
b=xlsread(inputfilename,'Geometry','A1');
h=xlsread(inputfilename,'Geometry','B1');


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

nel=size(nconn,1)

nnodes=size(ncoor,1);

KK=zeros(2*nnodes);

for i=1:nel
    disp(i)
    X= ncoor(nconn(i,:),1)';
    
    k=Beam2Noded_1Dk(X,b,h,E);
    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)';
            s1=[2*nconn(i,:)-1; 2*nconn(i,:)];
            sgdof=s1(:);
            [BM, SF, qpts]=Beam2Noded_1D_BM_SF(X, U(sgdof), b, h, E);
            qptsx(i,:)=qpts;
            %sss(i,:,:)=s';
            BMG(i,:)=BM;
            SFG(i,:)=SF;
        end        
    end

    figure(1)
    clf;
    hold off
    title('Deformed shape and indication w and slope');
    xdata=reshape( ncoor(nconn(:,:),1), size(nconn,1),size(nconn,2));
    udata=reshape( U(index_ugdof(2*(nconn(:,:)-1)+1,1)), size(nconn,1),size(nconn,2));
    plot(xdata,udata, 'b')
    hold on
    vdata=reshape( U(index_ugdof(2*(nconn(:,:)-1)+2,1)), size(nconn,1),size(nconn,2));
    plot(xdata,vdata, 'g')
    
    figure(2)
    clf;
    hold off
    title('BM and SF');
    plot( reshape(qptsx(:,:)',size(qptsx(:),1),1), reshape(BMG(:,:)',size(BMG(:),1),1),'r')
    hold on
    plot( reshape(qptsx(:,:)',size(qptsx(:),1),1), reshape(SFG(:,:)',size(SFG(:),1),1),'y')
    hold off
    
    flagfv=0;
    flagbc=0;
    flagscale=0;
    
end


