clear all;
uu=1;
E=250000;
A=2.0108;
B=0.9714;
C=9.1412;
D=0.2312;
fc=250;
clc;

% Reading number of nodes,elements and materials

    numnodes=input('number of nodes ? ');
    
    numelements=input('number of elements ? ');
      

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

% Reading Materials

    clc;
    display(['         materials :']);
    fc=input('f(c) ? ');
    E=input('E = ? ');
    h=input('y ? ');
    I=input('I ? ');
    Area=input('A ? ');


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

numeqk=numnodes*2-numconstrains;

for i=1:numnodes
    for j=1:2
        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)]);
        iele(i,1)=input('Node(i) ? ');
        iele(i,2)=input('Node(j) ? ');
        iele(i,3)=1;
end

% Calculate & Assemble stiffness matrix

k=zeros(numeq+1);
for i=1:numelements
    L=abs(coordinate(iele(i,2),1)-coordinate(iele(i,1),1));
    U=Area/L;
    V=12*E*I/L^3;
    W=6*E*I/L^2;
    X=4*E*I/L;
    Z=2*E*I/L;
    ke=[V,W,-V,W;W,X,-W,X/2;-V,-W,V,-W;W,X/2,-W,X];
    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?');
    L=abs(coordinate(iele(i,2),1)-coordinate(iele(i,1),1));
    flocal=[q*L/2;q*(L.^2)/12;q*L/2;-q*(L.^2)/12];
    F([iadress(iele(a,1),:),iadress(iele(a,2),:)],1)=flocal+F([iadress(iele(a,1),:),iadress(iele(a,2),:)],1);
    d=input('Is it another distributed force?(y/n)','s');
end


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

elementforces=zeros(numelements,6);
for i=1:numelements
    
    
    L=abs(coordinate(iele(i,2),1)-coordinate(iele(i,1),1));
    d1=U1(iadress(iele(i,1),1),1);
    d2=U1(iadress(iele(i,1),2),1);
    d3=U1(iadress(iele(i,2),1),1);
    d4=U1(iadress(iele(i,2),2),1);
    
    for x=-1:.1:1;
        M=(E*I/L.^2)*(6*x*d1+(3*x-1)*L*d2-6*x*d3+(3*x+1)*L*d4);
        V=(6*E*I/L.^3)*(2*d1+L*d2-2*d3+L*d4);
        siqma=M*h/I;
        tau=V/Area;
        J2=(siqma.^2+3*tau.^2)/3;
        siqma1=.5*siqma+(.25*siqma.^2+tau.^2).^.5;
        uu=uu+1;
        if i==1
            aqqq=L;
        end
        aq(uu,1)=(x+1)*L*.5+coordinate(iele(i,1),1);
        %aq(uu,1)=x;
        %aq(uu,2)=coordinate(iele(i,1),1)
        aq(uu,2)=M;
        aq(uu,3)=V;
        aq(uu,4)=siqma;
        aq(uu,5)=tau;
        aq(uu,6)=J2;
        aq(uu,7)=(Area*J2/fc.^2)+B*J2.^.5/fc+C*siqma1+D*siqma/fc-1;
        aq(uu,2)=M;
        aq(uu,3)=V;
    end
    
end