% ---------- BEAM ---------------
% Beam example solved using FEM
%   working on octave 3.6.2
% -------------------------------
%
% License:
%   Distributed under Creative Commons License
%   Attribution-ShareAlike 3.0 Unported
%
% Description:
%   In this example is solved the 
%   problem of a beam with zero
%   displacement and rotation of the
%   left node and loads applied at the 
%   right node. One beam element is used
%   with third degree shape functions.
%   Loads are given by the parameters P and M
%
% Updates:
%   2012-set : creation
%   2012-dic : several comments added
%
% -------------------------------
% authors:
%   J.M. Perez Zerpa - www.fing.edu.uy/~jperez/ 
%   P Castrillo      - www.fing.edu.uy/~pabloc/ 
%   IET, UdelaR, Set 2012
%
% -------------------------------
% -------------------------------

function fem_beam_example

close all,

% --- parameters ---
l = 2 ;          % length

P = -2 ;          % upward vertical load at right node
M = 1 ;          % counter-clockwise moment at right node
E = 1 ;          % Young modulus
I = 1 ;          % Second order Inertia of the beam 


% here the nodal forces are imposed
F = [  0 0 P  M ]' ;


% --- calculations... ---
%
% sitffness matrix
Ke = E*I/l^3 * [ ...
  12   6*l    -12   6*l   ;
  6*l  4*l^2  -6*l  2*l^2 ;
  -12  -6*l    12   -6*l   ;
  6*l  2*l^2  -6*l  4*l^2  ]
%
K = Ke ;
%
gdln         = [ 1 2 ]' ;
K(gdln , : ) = [] ;
K(:, gdln)   = [] ;
F(gdln)      = [] ;

q = K \ F 

REACC = Ke([1 2],[3 4]) * q

q = [ 0 0 q(1) q(2) ]' ;



% ----------------------------------------------
% --------  plotting  --------------------------
%
pasox = 0.1 ;   % step for graphics of solicitations
%
% vector of posotions for generating plots of solicitations
x     = 0 : pasox : l ;
O     = zeros(size(x)); % auxiliar vector for plots


% --- vertical displacement plotting -----------
%
v = zeros(size(x)); % vertical displacements

for i=1:length(v),
  %
  N = forma(x(i),l);
  v(i) = N * q ;
end

plot(x,v,'b-x')
hold on, axis equal
grid on
plot(x,O,'k-o')
title('Vertical displacement')
%~ axis([-l*.1 l*1.1 min(min(O),min(v)) max(max(O),max(v))] )
print('vertic.png','-dpng')

% --- counterclockwise rotation plotting -----------
%
t = zeros(size(x)); % angle of rotation

for i=1:length(t),
  %
  N = formapri(x(i),l);
  t(i) = N * q ;
end

figure
plot(x,t,'b-x')
grid on
hold on
plot(x,O,'k-o')
title('Rotation')
%~ axis([-l*.1 l*1.1 min(min(O),min(t)) max(max(O),max(t))] )
print('rotat.png','-dpng')




% --- momento
m = zeros(size(x));

for i=1:length(m),
  %
  N = formaseg(x(i),l);
  m(i) = N * q * E*I ;
end

figure
plot(x,m,'b-x')
grid on
hold on
plot(x,O,'k-o')
title('Moments')
%~ axis([-l*.1 l*1.1 min(min(O),min(m)) max(max(O),max(m))] )
print('momen.png','-dpng')





% --- Shear ---
V = zeros(size(x));

for i=1:length(V),
  %
  N = formater(x(i),l);
  V(i) = N * q * E*I ;
end

figure
plot(x,V,'b-x')
grid on
hold on
plot(x,O,'k-o')
title('Shear')
%~ axis([-l*.1 l*1.1 min(min(O),min(V)) max(max(O),max(V))] )
print('shear.png','-dpng')

% --- shape functions ---
function N = forma (x , l)

N1 =  (2*x^3 - 3*l*x^2 + l^3 )/ l^3 ;
N2 =  (x^3 -2 * l * x^2 + l^2*x  )/ l^2 ;
N3 = -( 2*x^3 - 3*l*x^2 )/l^3 ;
N4 =  (x^3 - l*x^2 )/l^2 ;

N= [N1 N2 N3 N4] ;


% --- shape functions first derivative ---
function N = formapri(x , l)

N1 =  (6*x^2 - 6*x*l )  / l^3 ;
N2 =  (3*x^2 -4*l*x +l^2)   / l^2 ;
N3 = -(6*x^2 - 6*x*l) /l^3 ;
N4 =  (3*x^2 - 2*l*x ) /l^2 ;

N= [N1 N2 N3 N4] ;


% --- shape functions second derivative ---
function N = formaseg (x , l)

N1 =  (12*x-6*l )  / l^3 ;
N2 =  (6*x-4*l )   / l^2 ;
N3 = -(12*x - 6*l) /l^3 ;
N4 =  (6*x - 2*l ) /l^2 ;

N= [N1 N2 N3 N4] ;


% --- shape functions third derivative ---
function N = formater (x , l)

N1 =  (12)  / l^3 ;
N2 =  (6 )   / l^2 ;
N3 = -(12) /l^3 ;
N4 =  (6 ) /l^2 ;

N= [N1 N2 N3 N4] ;

