%% Beam2D: Deflection of a beam using different types of beam elements
% The deflection of a beam of the length l, fixed on one side and charged with a load
% Q on the free side is calculated using different kinds of beam elements. (Lyonel
% Ehrl, 2006, for the seminar 'The Finite Element Method and the Analysis of Systems
% with Uncertain Properties')
close all; clear all; clc;

%% Independent Parameters
% Geometric parameters
    % Length l / m
    l = 0.28;
    % Height h / m
    h = .06;
    % Thickness b / m
    b = .004;
    % Plotting factor / -
    f = 100;
%%    
% Material properties
    % Young modulus E / Pa
    E = 2e11;
    % Poisson's ratio nu / -
    nu = 0.3;
%%
% Boundary conditions
    % Load at length l / N
    Q = -1e3;
    % Displacement at length 0 / m
    u = 0;
    % Deriverative of the displacement at length 0 / m/m
    du = 0;
%%    
% Simulation parameters
    % Shear correction factor k / -
    k = 5/6;
    % Number of elements nE / -
    nE = 10;
    % Beam element type / -
    beamtype = 4;
        % 1: 2-node Timoshenko
        % 2: 2-node Assumed Natural Strain (ANS)
        % 3: 2-node Shear deformation not considered
        % 4: 2-node exact Timoshenko
    ElementName = cell(1,4);
    ElementName = {' Timoshenko',' Assumed Natural Strain',' Only Bending',' exact Timoshenko'};
  
%% Dependent Parameters
% Geometrical moment of inertia / m4
I = b*h^3/12;
%%
% Area of beam crosssection A / m2
A = b*h;
%%
% Shear modulus G / Pa
G = E/2/(1+nu);
%%
% Parameters for exact Timoshenko shape functions
AL = k*G*A*(l/nE)^2/(E*I);
BE = 12 + AL;
%%
% Analytical model including shear deformation to calculate displacement UA / cm
UA = (1/3*Q*l^3/(E*I)+Q*l/(5/6*A*G))*100;
%%
% Basic stiffness matrix entries
    % for bending EB / -
    if (beamtype==1)||(beamtype==2),
        EB = [ 0  0  0  0 ;
               0  1  0 -1 ;
               0  0  0  0 ;
               0 -1  0 +1 ];
    elseif (beamtype==3),
        EB = [ 12        6*(l/nE)   -12         6*(l/nE) ;
               6*(l/nE)  4*(l/nE)^2  -6*(l/nE)  2*(l/nE)^2 ;
               -12      -6*(l/nE)    12        -6*(l/nE) ;
               6*(l/nE)  2*(l/nE)^2  -6*(l/nE)  4*(l/nE)^2 ];
    elseif (beamtype==4),
        SM = [  12/(l/nE)^3  6/(l/nE)^2        -12/(l/nE)^3  6/(l/nE)^2 ;
                 6/(l/nE)^2  4/(l/nE)*(1+3/AL)  -6/(l/nE)^2  2/(l/nE)*(1-6/AL) ;
               -12/(l/nE)^3 -6/(l/nE)^2         12/(l/nE)^3 -6/(l/nE)^2 ;
                 6/(l/nE)^2  2/(l/nE)*(1-6/AL)  -6/(l/nE)^2  4/(l/nE)*(1+3/AL) ];
    else 
        disp('Wrong beam type');
    end
    % for shearing ES / -
    if beamtype==1,
        ES = [    1/(l/nE) -1/2          -1/(l/nE) -1/2 ;
               -1/2         1/3*(l/nE)  1/2         1/6*(l/nE) ;
                 -1/(l/nE)  1/2           1/(l/nE)  1/2 ;
               -1/2         1/6*(l/nE)  1/2         1/3*(l/nE) ];
    elseif beamtype==2,
        ES = [    1/(l/nE) -1/2          -1/(l/nE) -1/2 ;
               -1/2         1/4*(l/nE)  1/2         1/4*(l/nE) ;
                 -1/(l/nE)  1/2           1/(l/nE)  1/2 ;
               -1/2         1/4*(l/nE)  1/2         1/4*(l/nE) ];
    elseif (beamtype==3)||(beamtype==4),
        ES = zeros(4);
    else 
        disp('Wrong beam type');
    end
%%    
% Prefactor for stiffness matrix entries
    % for bending PB / Nm
    if (beamtype==1)||(beamtype==2),
        PB = E*I/(l/nE);
    elseif (beamtype==3),
        PB = E*I/(l/nE)^3;
    elseif (beamtype==4),
        PB = 0;
    else 
        disp('Wrong beam type');
    end
    % for shearing PS / N
    PS = G*A*k;

%% Element Stiffness Matrix KE
    % Initializing element stiffness matrix KE
    KE = zeros(4,4);
    % Applying beam element type on element stiffness matrix KE
    if (beamtype==1)||(beamtype==2)||(beamtype==3),
        KE = PB.*EB + PS.*ES;
    elseif (beamtype==4),
        KE = AL*(E*I)/BE*SM;
    else 
        disp('Wrong beam type');
    end
    
%% Global Stiffness Matrix KG
% Initializing global stiffness matrix KG
KG = sparse(zeros(4+(nE-1)*2));
%%
% Applying element stiffness matrices KE on global stiffness matrix KG
for i=1:nE,
    r = 1+(i-1)*2;
    KG(r:r+3,r:r+3) = KG(r:r+3,r:r+3) + KE(1:4,1:4);
end
   
%% Calculating inverse global stiffness matrix KGinv (only necessary part)
KGinv = inv(KG(3:end,3:end));

%% Global load vector RG / N
if (beamtype==1)||(beamtype==2)||(beamtype==3)||(beamtype==4),
    RG = Q*[zeros((nE-1)*2,1);1;0];
else 
    disp('Wrong beam type');
end

%% Global Displacement Vector UG / m
UG = KGinv*RG;

%% Displacement at Length l, Ul / cm
Ul = UG(end-1:end)*100;

%% Plot Results
    % Plot Deformation
    figure(1);
    set(1,'Position',get(0,'ScreenSize'));
    x = linspace(0,l,(nE+1))';
    y = UG(1:2:end);
    dydx = diff([0;y])./diff(x);
    alph = atan(dydx);
    if (beamtype==1)||(beamtype==2)||(beamtype==4),
        beta = UG(2:2:end);
    elseif (beamtype==3),
        beta = UG(2:2:end)*0;
    end
    iniBeamx = [0 0 l l 0];
    iniBeamy = [-h/2 h/2 h/2 -h/2 -h/2];
    defBeamx = [0;0;(x-h/2*[0;sin(alph+beta)]*f);(flipud(x)+h/2*[sin(alph+beta);0]*f)];
    defBeamy = [-h/2;h/2;h/2;...
                (y*f+h/2+sign(alph+beta).*(h/2*(1-cos(alph+beta))*f));...
                (flipud(y)*f-h/2+sign(alph+beta).*(h/2*(1-cos(alph+beta))*f));-h/2];
    maxx = max(defBeamx)*1.2;
    miny = min(defBeamy)*1.2;
    maxy = max(defBeamy)*1.2;

    plot(iniBeamx,iniBeamy,'r-',...
         defBeamx,defBeamy,'k-o');
    set(findobj('Type','Line'),'LineWidth',3.0);
    set(gca,'FontSize',28,'XGrid','on','YGrid','on');
    title(strcat('Deformation Profile, ',ElementName{beamtype},', ',num2str(nE),' elements'),'FontSize',28);
    xlabel('x / m','FontSize',28);
    ylabel('y / m','FontSize',28);
    axis([0 maxx -maxx maxy ]);
    legend('undeformed',strcat('deformed*',num2str(f)),3);