% Introduction of Finite Elements in Engineering Mechanics 
% Dr. X. Frank Xu, Stevens Institute of Technology

clear
%----------------------
%Mesh Information
%----------------------
NumNP = 4; %number of nodal points
NumEL = 3; % number of elements
NP=[1 2;2 3;3 4]; %Nodes in elements from 1 to NumEL
XORD(1:NumNP,1)=[0,240,240,240]'; % X-coordinate of nodal points
YORD(1:NumNP,1)=[0,0,72,144]'; % Y-coordinate of nodal points

%----------------------
%Specify Boundary Conditions
%----------------------
NumEQ = 3*NumNP;
U=zeros(NumEQ,1);
F=zeros(NumEQ,1);
NPBC=zeros(NumEQ,1);
NPBC(1)=1; %X-Displacement of Node 1 is constrained
NPBC(2)=1; 
NPBC(3)=1; 
NPBC(10)=1; 
NPBC(11)=1; 
NPBC(12)=1; 
U(1)=0; % Specified X-displacement for Node 1 
U(2)=0;  
U(3)=0;  
U(10)=0;  
U(11)=0; 
U(12)=0;  

F(5)=-3*20/2; 
F(6)=3*240^2/12/12; 
F(7)=-10; 

%----------------------
%Specify Material Properties
%----------------------
E(1:NumEL,1)=29000; % Young's modulus (unit:ksi)
A(1:NumEL,1)=[12;8;8]; % Cross section area(unit:in2) 
II(1:NumEL,1)=[600;300;300]; % Moment Inertia(unit:in4) 
%----------------------
% Form Stiffness Matrix
%----------------------
K=zeros(NumEQ, NumEQ); 
for I=1:NumEL
    Ke=zeros(6,6);
    L(I)=((XORD(NP(I,1))-XORD(NP(I,2)))^2+(YORD(NP(I,1))-YORD(NP(I,2)))^2)^0.5;
    BarAngle(I)=atan2((YORD(NP(I,2))-YORD(NP(I,1))),(XORD(NP(I,2))-XORD(NP(I,1))));
      C=cos(BarAngle(I));
    S=sin(BarAngle(I));
    Ke(1,:)=E(I)/L(I)*[A(I)*C*C+12*II(I)*S^2/L(I)^2  C*S*(A(I)-12*II(I)/L(I)^2)    -6*II(I)*S/L(I) -A(I)*C*C-12*II(I)*S*S/L(I)^2 (-A(I)+12*II(I)/L(I)^2)*S*C    -6*II(I)*S/L(I)];
    Ke(2,:)=E(I)/L(I)*[C*S*(A(I)-12*II(I)/L(I)^2)    A(I)*S*S+12*II(I)*C^2/L(I)^2  6*II(I)*C/L(I)  -A(I)*S*C+12*II(I)*C*S/L(I)^2 -A(I)*S*S-12*II(I)*C*C/L(I)^2  6*II(I)*C/L(I)];
    Ke(3,:)=E(I)/L(I)*[-6*II(I)*S/L(I)               6*II(I)*C/L(I)                4*II(I)         6*II(I)*S/L(I)                -6*II(I)*C/L(I)                2*II(I)];
    Ke(4,:)=E(I)/L(I)*[-A(I)*C*C-12*II(I)*S*S/L(I)^2 -A(I)*S*C+12*II(I)*C*S/L(I)^2 6*II(I)*S/L(I)  A(I)*C*C+12*II(I)*S^2/L(I)^2  A(I)*C*S-12*II(I)*S*C/L(I)^2   6*II(I)*S/L(I)];
    Ke(5,:)=E(I)/L(I)*[-A(I)*C*S+12*II(I)*S*C/L(I)^2 -A(I)*S*S-12*II(I)*C*C/L(I)^2 -6*II(I)*C/L(I) A(I)*C*S-12*II(I)*S*C/L(I)^2  A(I)*S*S+12*II(I)*C*C/L(I)^2   -6*II(I)*C/L(I)];
    Ke(6,:)=E(I)/L(I)*[-6*II(I)*S/L(I)               6*II(I)*C/L(I)                 2*II(I)        6*II(I)*S/L(I)                -6*II(I)*C/L(I)                 4*II(I)]; %%Element stiffness matrix
    for J=1:2  
        PosK(J*3-2)=NP(I,J)*3-2;% Position of Ke in global stiffness matrix K 
        PosK(J*3-1)=NP(I,J)*3-1;
        PosK(3*J)=NP(I,J)*3; % Position of Ke in global stiffness matrix K 
    end
    K(PosK,PosK)=K(PosK,PosK)+Ke; % Global matrix 
end
K0=K;
%--------------------------------------
% IMPLEMENTATION OF BOUNDARY CONDITIONS
%--------------------------------------
for I=1:NumEQ
    if NPBC(I)==1 
       K(I,I)=K(I,I)*1E9; % Big number approach
       F(I)=K(I,I)*U(I);
    end
end

%----------------------
% Solve displacement
%----------------------
U = K\F;

%-----------------------------------------
% Post-processing to obtain reaction forces
%-----------------------------------------

FBC=K0*U;
       
%---------------------------------------------
% Plot initial shape of the beam
ux0=XORD; vy0=YORD;
plot(ux0,vy0,'-o'); hold on

% -------------------------------------------------------------------------
% Plot the enlarged deformation of the beam elements
%
 r=300; % deformation enlargement ratio
 N=30;
for i=1:NumEL
    [xord,yord,x1,y1,R,ux_new,vy_new] = frame_plot(U,i,XORD,YORD,NumNP,NP,r,N);
end
title('Deformation of the frame (ratio=300)')