clear all
close all

% ------ Generate the model ------------------------------------------
% ------ material data ------------------------------------------
E=2.5e10; rho=2400;
b=0.6
h=0.8
A = b*h
I = (b*h^3)/12
ep=[E A I rho*A]
% ------ topology -----------------------------------------------
Edof=[1   1  2  3  4  5  6
      2   4  5  6  7  8  9
      3   7  8  10 11 12 13
      4   11 12 13 14 15 16
      5   14 15 16 17 18 19
      6   17 18 19 20 21 22]
% ------ list of coordinates  -----------------------------------
Ex=[0 3
    3 6
    6 10
    10 14
    14 14
    14 14]

Ey=[0 4
    4 8
    8 8
    8 8
    8 4
    4 0]

% ----------------------------------------------------------------------
% Opole: 14-12-2013 17:34:17
%
% ----------------------------------------------------------------------
[n_el,xxx] = size(Edof);
dofmax = max(max(Edof));

% ------ generate element matrices, assemble in global matrices -
K=zeros(dofmax);     M=zeros(dofmax);
for i=1:n_el
  [k,m,c]=beam2d(Ex(i,:),Ey(i,:),ep);
  K=assem(Edof(i,:),K,k);  M=assem(Edof(i,:),M,m);
end

% ----- Draw a plot of the element mesh --------------------------
clf;     eldraw2(Ex,Ey,[1 2 2],Edof);
grid;    title('Geometria ukladu nieodksztalconego')

% ----- Eigenvalue analysis --------------------------------------
% - Warunki brzegowe -
b=[1 2 3 20 21 22]';

% ----------------------------------------------------------------------
% Opole: 20-04-2013 15:12:44
% metoda mas skupionych
% ----------------------------------------------------------------------
% M = eye(size(K))
% M(1,1)=1
% M(2,2)=1
% M(4,4)=5760
% M(5,5)=5760
% M(7,7)=5184
% M(8,8)=5184
% M(11,11)=4608
% M(12,12)=4608
% M(14,14)=4608
% M(15,15)=4608
% M(17,17)=4608
% M(18,18)=4608
% M(20,20)=1
% M(21,21)=1

[La,Egv]=eigen(K,M,b);
Omega=sqrt(La); % rad/s
Freq=Omega/(2*pi); % Hz
Okres=1./Freq; % s

% ----- plot first eigenmode ---------------------------------------

figure,    clf,     grid,
title(['Pierwsza postac: F_1 = ',num2str(Freq(1)),'Hz, \omega_1 = ' ...
       ,num2str(Omega(1)),'rad/s, T_1 = ',num2str(Okres(1)), 's']),
eldraw2(Ex,Ey,[2 3 1]);
Edb=extract(Edof,Egv(:,1));      eldisp2(Ex,Ey,Edb,[1 2 2]);

% ----- plot second eigenmode ---------------------------------------

figure,    clf,     grid,
title(['Druga postac: F_2 = ',num2str(Freq(2)),'Hz, \omega_2 = ' ...
       ,num2str(Omega(2)),'rad/s, T_2 = ',num2str(Okres(2)), 's']),
eldraw2(Ex,Ey,[2 3 1]);
Edb=extract(Edof,Egv(:,2));      eldisp2(Ex,Ey,Edb,[1 2 2]);

% ----- plot third eigenmode ---------------------------------------

figure,    clf,     grid,
title(['Trzecia postac: F_3 = ',num2str(Freq(3)),'Hz, \omega_3 = ' ...
       ,num2str(Omega(3)),'rad/s, T_3 = ',num2str(Okres(3)),'s']),
eldraw2(Ex,Ey,[2 3 1]);
Edb=extract(Edof,Egv(:,3));      eldisp2(Ex,Ey,Edb,[1 2 2]);

% ----- plot fourth eigenmode ---------------------------------------

figure,    clf,     grid,
title(['Czwarta postac: F_4 = ',num2str(Freq(4)),'Hz, \omega_4 = ' ...
       ,num2str(Omega(4)),'rad/s, T_4 = ',num2str(Okres(4)),'s']),
eldraw2(Ex,Ey,[2 3 1]);
Edb=extract(Edof,Egv(:,4));      eldisp2(Ex,Ey,Edb,[1 2 2]);

% ----- plot fifth eigenmode ---------------------------------------

figure,    clf,     grid,
title(['Piata postac: F_5 = ',num2str(Freq(5)),'Hz, \omega_5 = ' ...
       ,num2str(Omega(5)),'rad/s, T_5 = ',num2str(Okres(5)),'s']),
eldraw2(Ex,Ey,[2 3 1]);
Edb=extract(Edof,Egv(:,5));      eldisp2(Ex,Ey,Edb,[1 2 2]);

% ----- plot eight eigenmodes ------------------------------------

% figure, clf, axis('equal'), hold on, axis off
% sfac=0.5;
% title('The first eight eigenmodes (Hz)' )
% for i=1:4;
%   Edb=extract(Edof,Egv(:,i));
%   Ext=Ex+(i-1)*3;                eldraw2(Ext,Ey,[2 3 1]);
%   eldisp2(Ext,Ey,Edb,[1 2 2],sfac);
%   FreqText=num2str(Freq(i));     text(3*(i-1)+.5,1.5,FreqText);
% end;
% Eyt=Ey-4;
% for i=5:8;
%   Edb=extract(Edof,Egv(:,i));
%   Ext=Ex+(i-5)*3;                eldraw2(Ext,Eyt,[2 3 1]);
%   eldisp2(Ext,Eyt,Edb,[1 2 2],sfac);
%   FreqText=num2str(Freq(i));     text(3*(i-5)+.5,-2.5,FreqText);
% end

% -------------------- end ---------------------------------------
echo off
