% 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','A:A');
h=xlsread(inputfilename,'Geometry','B:B');


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)';
    Y= ncoor(nconn(i,:),2)';
    k=planetruss2Noded(X,Y,b(i),h(i),E);
    s1=[2*nconn(i,:)-1; 2*nconn(i,:)];
    sgdof=s1(:);
    
    KK(sgdof, sgdof)=KK(sgdof,sgdof) + k;
end
    F=zeros(nnodes*2,1);

    deformscaleX=xlsread(inputfilename,'scale','A1');
    deformscaleY=xlsread(inputfilename,'scale','B1');
    
    fv=xlsread(inputfilename,'force');
    nfv=size(fv,1);
    F((fv(:,1)-1)*2+fv(:,2))=fv(:,3) 
    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)
    
    
    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)';
        L= (  (X(end)-X(1))^2 + (Y(end)-Y(1))^2 )^0.5;
        s1=[2*nconn(i,:)-1; 2*nconn(i,:)];
        X1=X + U(s1(1,:))';
        Y1=Y + U(s1(2,:))';
        L1= (  (X1(end)-X1(1))^2 + (Y1(end)-Y1(1))^2 )^0.5;
        delta=L1-L;
        P=b(i).*h(i)*E*delta/L;
        s=E*delta/L;
        PG(i,:)=P;
        SG(i,:)=s;
    end        
    
    figure(1)
    clf;
    hold off
    title('Deformed shape and indication of member forces');
    xdata=reshape( ncoor(nconn(:,:),1), size(nconn,1),size(nconn,2));
    ydata=reshape( ncoor(nconn(:,:),2), size(nconn,1),size(nconn,2));
    udata=reshape( U(index_ugdof(2*(nconn(:,:)-1)+1,1)), size(nconn,1),size(nconn,2));
    vdata=reshape( U(index_ugdof(2*(nconn(:,:)-1)+2,1)), size(nconn,1),size(nconn,2));
    p1=patch(xdata', ydata', (ydata>-1)'*2)
    set(p1,'EdgeColor','b')
    set(p1,'EdgeColor','flat', 'LineWidth',3)
    hold on
    p2=patch((xdata+udata)', (ydata+vdata)', (ydata>-1)'*2)
    set(p2,'EdgeColor','r')
    Htxt=text([sum(xdata(:,:),2)/2-.05], [sum(ydata(:,:),2)/2+.05], [num2str(PG)], 'HorizontalAlignment','center')
    figure(2)
    clf;
    hold off
    title('Deformed shape and indication of stress in members');
    xdata=reshape( ncoor(nconn(:,:),1), size(nconn,1),size(nconn,2));
    ydata=reshape( ncoor(nconn(:,:),2), size(nconn,1),size(nconn,2));
    udata=reshape( U(index_ugdof(2*(nconn(:,:)-1)+1,1)), size(nconn,1),size(nconn,2));
    vdata=reshape( U(index_ugdof(2*(nconn(:,:)-1)+2,1)), size(nconn,1),size(nconn,2));
    p1=patch(xdata', ydata', (ydata>-1)'*2)
    set(p1,'EdgeColor','b')
    set(p1,'EdgeColor','flat', 'LineWidth',3)
    hold on
    p2=patch((xdata+udata)', (ydata+vdata)', (ydata>-1)'*2)
    set(p2,'EdgeColor','r')
    Htxt=text([sum(xdata(:,:),2)/2-.05], [sum(ydata(:,:),2)/2+.05], [num2str(SG)], 'HorizontalAlignment','center')