clc;
clear all;

% Reading number of nodes,elements and materials

while 1==1
    numnodes=input('number of nodes ? ');
    if numnodes<=1
        display('number of nodes must be more than 1 node');
    else break;
    end
end
while 1==1
    numelements=input('number of elements ? ');
    if numelements<numnodes-1
        display('number of elements must be more than number of nodes-2');
    else break;
    end
end
while 1==1
    nummaterials=input('number of materials ? ');
    if nummaterials<1
        display('number of materials must be more than 0 material');
    else break;
    end
end

coordinate=zeros(numnodes,2);
propma=zeros(nummaterials,4);
iele=zeros(numelements,3);
iadress=zeros(numnodes,3);
numeq=0;
numconstrains=0;
Res=zeros(numnodes,3);
Restrains=zeros(numnodes,3);

% Reading Materials
for i=1:nummaterials
    clc;
    display(['         material :',num2str(i)]);
    propma(i,1)=input('E ? ');
    propma(i,2)=input('I ? ');
    propma(i,3)=input('A ? ');
end

% Reading coordinates & Boundry condition of nodes
for i=1:numnodes
    clc;
    fprintf('X coordinate of node %1.0f',i);
    coordinate(i,1)=input('? ');
    fprintf('Y coordinate of node %1.0f',i);
    coordinate(i,2)=input('? ');
    while 1==1
        Res(i,1)=input('X Restrain(0 for free,1 for closed) ? ');
        if Res(i,1)>=2
            display('X Restrain must be 0 for free,1 for closed');
        else break;
        end
    end
    while 1==1
        Res(i,2)=input('Y Restrain(0 for free,1 for closed) ? ');
        if Res(i,2)>=2
            display('Y Restrain must be 0 for free,1 for closed');
        else break;
        end
    end
    while 1==1
        Res(i,3)=input('rotation Restrain(0 for free,1 for closed) ? ');
        if Res(i,3)>=2
            display('rotation Restrain must be 0 for free,1 for closed');
         else break;
        end
    end
    numconstrains=numconstrains+Res(i,1)+Res(i,2)+Res(i,3);
    % iadress(i,:)=[3*i-2,3*i-1,3*i];
    
end

numeqk=numnodes*3-numconstrains;

for i=1:numnodes
    for j=1:3
        if Res(i,j)==0
            numeq=numeq+1;
            iadress(i,j)=numeq;
        else
            iadress(i,j)=numeqk+1;
        end
    end
end

% Reading connectivity of elements

for i=1:numelements
    clc;
    display(['      Element : ',num2str(i)]);
    while 1==1
        iele(i,1)=input('Node(i) ? ');
        if iele(i,1)<=0 || iele(i,1)>numnodes
            display('This node does not exist,input correct one again');
        else break;
        end
    end
    while 1==1
        iele(i,2)=input('Node(j) ? ');
        if iele(i,2)<=0 || iele(i,2)>numnodes
            display('This node does not exist,input correct one again');
        else break;
        end
    end
    while 1==1
        iele(i,3)=input('Material of element ? ');
        if iele(i,3)<=0 || iele(i,3)>nummaterials
            display('This materal does not exist,input correct one again');
        else break;
        end
    end    
end

% Calculate & Assemble stiffness matrix

k=zeros(numeq+1);
for i=1:numelements
    deltax=coordinate(iele(i,2),1)-coordinate(iele(i,1),1);
    deltay=coordinate(iele(i,2),2)-coordinate(iele(i,1),2);
    L=(deltax^2+deltay^2)^.5;
    cx=deltax/L;
    sx=deltay/L;
    E=propma(iele(i,3),1);
    I=propma(iele(i,3),2);
    A=propma(iele(i,3),3);
    U=E*A/L;
    V=12*E*I/L^3;
    W=6*E*I/L^2;
    X=4*E*I/L;
    Z=2*E*I/L;
    R=[cx sx 0 0 0 0;-sx cx 0 0 0 0;0 0 1 0 0 0;0 0 0 cx sx 0;0 0 0 -sx cx 0;0 0 0 0 0 1];
    ke=[U 0 0 -U 0 0;0 V W 0 -V W;0 W X 0 -W Z;-U 0 0 U 0 0;0 -V -W 0 V -W;0 W Z 0 -W X];
    Ke=R'*ke*R;
    nodei=iele(i,1);
    nodej=iele(i,2);
    k([iadress(nodei,:) iadress(nodej,:)],[iadress(nodei,:),iadress(nodej...
        ,:)])=k([iadress(nodei,:) iadress(nodej,:)],[iadress(nodei,:),...
        iadress(nodej,:)])+Ke;
end

% Reading concentrated forces and calculate force matrix

clc;
F=zeros(numeq+1,1);
c=input('is it any concentrated force or moment?(y/n)','s');
while c=='y'
    a=input('Which node?');
    b=input('specify direction(1 for X,2 for Y,3 for moment)=');
    F(iadress(a,b),1)=input('How much?');
    c=input('Is it another concentrated force or moment?(y/n)','s');
end

% Reading distributed forces and calculate force matrix 

d=input('Is it any distributed force?(y/n)','s');
while d=='y'
    a=input('Which element?');
    q=input('How much?');
    deltax=coordinate(iele(a,2),1)-coordinate(iele(a,1),1);
    deltay=coordinate(iele(a,2),2)-coordinate(iele(a,1),2);
    cx=deltax/L;
    sx=deltay/L;
    L=(deltax^2+deltay^2)^.5;
    R=[cx sx 0 0 0 0;-sx cx 0 0 0 0;0 0 1 0 0 0;0 0 0 cx sx 0;0 0 0 -sx cx 0;0 0 0 0 0 1];
    flocal=[0;q*L/2;q*(L.^2)/12;0;q*L/2;-q*(L.^2)/12];
    
   Restrains(iele(a,1),:)=-flocal(1:3,1)'.*Res(iele(a,1),:)+Restrains(iele(a,1),:);
   Restrains(iele(a,2),:)=-flocal(4:6,1)'.*Res(iele(a,2),:)+Restrains(iele(a,2),:);
    
    
    %F2([iadress(iele(a,1),:) iadress(iele(a,2),:)],1)=R'*flocal;
    F([iadress(iele(a,1),:),iadress(iele(a,2),:)],1)=flocal+ F([iadress(iele(a,1),:),iadress(iele(a,2),:)],1);
    %F(iadress(iele(a,1:2),:),1)=flocal+F(iadress(iele(a,1:2),:),1)
    d=input('Is it another distributed force?(y/n)','s');
end

% calculate total force matrix

%F=F1+F2;

% calculate displacement matrix

U1=k(1:numeq,1:numeq)\F(1:numeq,1);
U1(numeq+1,1)=0;

elementforces=zeros(numelements,6);
for i=1:numelements
    
    deltax=coordinate(iele(i,2),1)-coordinate(iele(i,1),1);
    deltay=coordinate(iele(i,2),2)-coordinate(iele(i,1),2);
    L=(deltax^2+deltay^2)^.5;
    cx=deltax/L;
    sx=deltay/L;
    E=propma(iele(i,3),1);
    I=propma(iele(i,3),2);
    A=propma(iele(i,3),3);
    U=E*A/L;
    V=12*E*I/L^3;
    W=6*E*I/L^2;
    X=4*E*I/L;
    Z=2*E*I/L;
    R=[cx sx 0 0 0 0;-sx cx 0 0 0 0;0 0 1 0 0 0;0 0 0 cx sx 0;0 0 0 -sx cx 0;0 0 0 0 0 1];
    ke=[U 0 0 -U 0 0;0 V W 0 -V W;0 W X 0 -W Z;-U 0 0 U 0 0;0 -V -W 0 V -W;0 W Z 0 -W X];
    Ke=R'*ke*R;
    elementforces(i,:)=(Ke*U1([iadress(iele(i,1),:),iadress(iele(i,2),:)],1))';
    % Calculating Restrains Forces
    
   Restrains(iele(i,1),:)=elementforces(i,1:3).*Res(iele(i,1),:)+Restrains(iele(i,1),:);
   Restrains(iele(i,2),:)=elementforces(i,4:6).*Res(iele(i,2),:)+Restrains(iele(i,2),:);
    
    
end