clc;
clear all;
numnodes=input('number of nodes ? ');
numelements=input('number of elements ? ');
nummaterials=input('number of materials ? ');

coordinate=zeros(numnodes,2);
k=zeros(3*numnodes);
propma=zeros(nummaterials,4);
iele=zeros(numelements,4);
localk=zeros(6);
iadress=zeros(numnodes,3);
numeq=0;
numconstrains=0;
R=zeros(numnodes,3);

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

% Reading coordinates & Boundry condition of nodes
for i=1:numnodes
    clc;
    fprintf('X coordinate of node %2.0f',i);
    coordinate(i,1)=input('? ');
    fprintf('Y coordinate of node %2.0f',i);
    coordinate(i,2)=input('? ');
    R(i,1)=input('X Restrain ? ');
    R(i,2)=input('X Restrain ? ');
    R(i,3)=input('X Restrain ? ');
    numconstrains=numconstrains+R(i,1)+R(i,2)+R(i,3);
    % iadress(i,:)=[3*i-2,3*i-1,3*i];
    
end
numeq=numnodes*3-numconstrains;
for i=1:numnodes
for j=1:3
    if R(i,j)==1
    iadress(i,j)=numeq+1;
else
    iadress(i,j)=3*i-3+j;
    end
end
end


% Reading connectivity of elements

for i=1:numelements
    clc;
    display(['      Element : ',num2str(i)]);
    iele(i,1)=i;
    iele(i,2)=input('Node(i) ? ');
    iele(i,3)=input('Node(j) ? ');
    iele(i,4)=input('Material of element ? ');
    
end

% Calculate & Assemble stiffness matrix

for i=1:numelements

    deltax=coordinate(iele(i,3),1)-coordinate(iele(i,2),1);
    deltay=coordinate(iele(i,3),2)-coordinate(iele(i,2),2);
    L=(deltax^2+deltay^2)^.5;
    cx=deltax/L;
    sx=deltay/L;
    E=propma(iele(i,4),2);
    A=propma(iele(i,4),3);
    I=propma(iele(i,4),4);
    U=E*A/L;
    V=12*E*I/L^3;
    W=6*E*I/L^2;
    localk=[U*cx^2+V*sx^2,U*sx*cx-V*sx*cx,-W*sx,-U*cx^2-V*sx^2,-U*cx*sx+...
        V*cx*sx,-W*sx;(U-V)*sx*cx,V*cx^2+U*sx^2,W*cx,(-U+V)*cx*sx,...
        -U*sx^2-V*cx^2,W*cx;-W*sx,W*cx,4*E*I/L,W*sx,-W*cx,2*E*I/L;...
        -U*cx^2-V*sx^2,(-U+V)*cx*sx,W*sx,U*cx^2+V*sx^2,(U-V)*cx*sx,...
        W*sx;(-U+V)*sx*cx,-U*sx^2-V*cx^2,-W*cx,(U-V)*sx*cx,V*cx^2+...
        U*sx^2,-W*cx;-W*sx,-W*cx,2*E*I/L,W*sx,-W*cx,4*E*I/L];
    nodei=iele(i,2);
    nodej=iele(i,3);
    k([iadress(nodei,:) iadress(nodej,:)],[iadress(nodei,:),iadress(nodej...
        ,:)])=k([iadress(nodei,:) iadress(nodej,:)],[iadress(nodei,:),...
        iadress(nodej,:)])+localk;
end