function [K_free, M_free, DOF] = FCN_beam3(a,b,L,E,ro,e,M_type)
%
% Funkcia sluzi na vypocet globalnej matice tuhosti a hmotnosti
% neuvazbeneho telesa pouzitim 2D nosnikoveho linearneho prvku s dvoma 
% uzlovymi bodmi. Kazdy uzlovy bod ma 3 stupne volnosti - x, y, phi.
%
% Inputs:
%
%   a   -   sirka nosnika [m]
%   b   -   vyska nosnika [m]
%   L   -   dlzka nosnika [m]
%   E   -   Modul pruznosti materialu v tahu [Pa]
%   ro  -   Hustota materialu [kg/m]
%   e   -   Pocet elementov [-]
%   M_type  -   Typ matice hmotnosti
%               M_type = 'cons' -   Matica konzistentnej hmotnosti
%               M_type = 'lump' -   Matica sustredenej hmotnosti
%
% Outputs:
%
%   K_free  -   Globalna matica tuhosti neuvazbeneho telesa
%   M_free  -   Globalna matica hmotnosti neuvazbeneho telesa
%   DOF     -   pocet stupnov volnosti diskretizovaneho neupevneneho telesa
%
% Example:
%
%   [K_free, M_free, DOF] = FCN_beam3(0.01,0.01,1,2.068e11,7830,3,'cons')
%
% Urcene pre vyuku predmetu MKP v dynamike struktur I, II
% Autor:    Vedelin Hok

% Real constant (Realne konstanty) ----------------------------------------
S = a*b;
Jz = (1/12)*a*b^3;

% Meshing (Konecnoprvkova siet) -------------------------------------------
% Pocet uzlov
    N = e + 1;
    % Pocet DOF jedneho uzla
    DOF_N = 3;
    % Pocet vsetkych DOF
    DOF = DOF_N*N;
    % Dlzka elementu  
    Le = L/e;

% Meshing processor +++++++++++++++++++++++++++++++++++++++++++++++++++++++

% MATICA TUHOSTI PRUTOVEHO PRVKU ------------------------------------------
Ke = [E*S/Le   0               0             -E*S/Le   0                0;
     0        12*E*Jz/(Le^3)   6*E*Jz/(Le^2)  0       -12*E*Jz/(Le^3)   6*E*Jz/(Le^2);
     0        6*E*Jz/(Le^2)    4*E*Jz/Le      0       -6*E*Jz/(Le^2)    2*E*Jz/Le;
    -E*S/Le   0                0              E*S/Le   0                0;
     0       -12*E*Jz/(Le^3)  -6*E*Jz/(Le^2)  0        12*E*Jz/(Le^3)  -6*E*Jz/(Le^2);
     0        6*E*Jz/(Le^2)    2*E*Jz/Le      0       -6*E*Jz/(Le^2)    4*E*Jz/Le];

% Zostavenie bunkoveho pola rozsirenych matic tuhosti
j = 1;
for i = 1:1:e
    Ke_ = zeros(DOF);
    Ke_(j:j+5,j:j+5) = Ke;
    j = j + 3;
    Ke_cell{i} = Ke_;
end

% Globalna matica tuhosti volneho telesa - 
% Sumacia jednotlivych buniek Ke_cell
K_free = Ke_cell{1};
for i = 1:1:(e-1)
    K_free = K_free + Ke_cell{i+1};
end

% MATICA HMOTNOSTI PRUTOVEHO PRVKU-----------------------------------------
m = strcmp(M_type,'lump');

% Matica sustredenej (lumped) hmotnosti prutoveho prvku
if m == 1   
    me = ro*S*Le;
    Ie = (1/12)*me*Le^2; iz = sqrt(Jz/S); Is = iz^2*me;
    
    Me = blkdiag(me/2, me/2,(Ie+Is)/2, me/2, me/2,(Ie+Is)/2);
     
% Matica konzistentnej (consistent) hmotnosti prutoveho prvku    
else        
    me = ro*S*Le;    % Hmotnost elementu
    c = me/420;
    Me = c*[140   0       0           70    0         0;
            0     156     22*Le       0     54       -13*Le;
            0     22*Le   4*(Le^2)    0     13*Le    -3*(Le^2);
            70    0       0           140   0         0;
            0     54      13*Le       0     156      -22*Le;
            0    -13*Le  -3*(Le^2)    0    -22*Le     4*(Le^2)];
end

% Zostavenie bunkoveho pola rozsirenych matic hmotnosti
j = 1;
for i = 1:1:e
    Me_ = zeros(DOF);
    Me_(j:j+5,j:j+5) = Me;
    j = j + 3;
    Me_cell{i} = Me_;
end

% Globalna matica hmotnosti volneho telesa - 
% Sumacia jednotlivych buniek Me_cell
M_free = Me_cell{1};
for i = 1:1:(e-1)
    M_free = M_free + Me_cell{i+1};
end
